What is wrong with my code? I have to find value of x(20 sec), y(20 sec) & z(20 sec). The correct answers are: x(20) = 0.704, y(20) = 0.665, z(20) = 0.246.

function a=ques()
D1=@(t,x,y,z) z-y/2; % dx/dt
D2=@(t,x,y,z) x/2-z/sqrt(2); %dy/dt
D3=@(t,x,y,z) y/sqrt(2)-x/2; %dz/dt
h=0.01;
t(1)=0;x(1)=1;y(1)=0;z(1)=0;
for i=1:1999
    t(i+1)=t(i)+h;
    k11=h * D1(t(i),      x(i),          y(i),         z(i));
    k21=h * D2(t(i),      x(i),          y(i),         z(i));
    k31=h * D3(t(i),      x(i),          y(i),         z(i));
    
    k12=h * D1(t(i) + h/2, x(i) + k11/2, y(i) + k21/2, z(i) + k31/2);
    k22=h * D2(t(i) + h/2, x(i) + k11/2, y(i) + k21/2, z(i) + k31/2);
    k32=h * D3(t(i) + h/2, x(i) + k11/2, y(i) + k21/2, z(i) + k31/2);
    
    k13=h * D1(t(i) + h/2, x(i) + k12/2, y(i) + k22/2, z(i) + k32/2);
    k23=h * D2(t(i) + h/2, x(i) + k12/2, y(i) + k22/2, z(i) + k32/2);
    k33=h * D3(t(i) + h/2, x(i) + k12/2, y(i) + k22/2, z(i) + k32/2);
    
    k14=h * D1(t(i) + h,   x(i) + k13,   y(i) + k23,   z(i) + k33);
    k24=h * D1(t(i) + h,   x(i) + k13,   y(i) + k23,   z(i) + k33);
    k34=h * D1(t(i) + h,   x(i) + k13,   y(i) + k23,   z(i) + k33);
    
    x(i+1)=x(i)+ (k11+2*k12+2*k13+k14)/6;
    y(i+1)=y(i)+ (k21+2*k22+2*k23+k24)/6;
    z(i+1)=z(i)+ (k31+2*k32+2*k33+k34)/6;
end
x(2000)
y(2000)
z(2000)

Answer 1

[darova]

mistake

[Dhrumil Patadia]

Oh..

but even then i don't get the correct answer

[darova]

Maybe your 'correct' answers are not correct

Check the solution using ode45

f = @(t,u) [u(3)-u(2)/2
    u(1)/2-u(3)/sqrt(2)
    u(2)/sqrt(2)-u(1)/2];
[t,u] = ode45(f,[0 20],[1 0 0]);
plot(t,u,'.')

Looks ok