Runge-Kutta 4th Order on a 2nd Order ODE

I have the following 2nd order ODE:

y('')=exp(2*t)*sin(t) - 2*y + 2*y(')

i used Runge-Kutta 4th Order to solve it with below script but I cannot get the correct output... is it the code or my calculus?

for exemple: the value of k11 shoud be -0.06 but MATLAB showed it as 0.0771

any one can help me?

f1 = @(t,y,z) z;

f2 = @(t,y,z) sin(t)*exp(2*t) - 2*y + 2*z;

t(1) = 0;

z(1) = -0.6;

y(1) = -0.4;

h = 0.1;

tfinal = 1;

m = (tfinal-t(1))/h;

for j=1:(m-1)

k11 = h*feval( f1 , t(j) , y(j) , z(j) );

k21 = h*feval( f2 , t(j) , y(j) , z(j) );

k12 = h*feval( f1 , t(j)+h/2 , y(j)+k11/2 , z(j)+k21/2);

k22 = h*feval( f2 , t(j)+h/2 , y(j)+k11/2 , z(j)+k21/2);

k13 = h*feval( f1 , t(j)+h/2 , y(j)+k12/2 , z(j)+k22/2);

k23 = h*feval( f2 , t(j)+h/2 , y(j)+k12/2 , z(j)+k22/2);

k14 = h*feval( f1 , t(j)+h , y(j)+k13 , z(j)+k23 );

k24 = h*feval( f2 , t(j)+h , y(j)+k13 , z(j)+k23 );

y(j+1) = y(j) + (k11 + 2*k12 + 2*k13 + k14)/6;

z(j+1) = z(j) + (k21 + 2*k22 + 2*k23 + k24)/6;

t(j+1) = t(j) + h;

end

disp(y);

disp(z);

Answer 1

[Torsten]

Comparing with a MATLAB integrator, the results look ok.

I don't know how you derived that "the value of k11 shoud be -0.06 but MATLAB showed it as 0.0771"

Maybe you used a different h.

f1     = @(t,y,z) z;
f2     = @(t,y,z) sin(t)*exp(2*t) - 2*y + 2*z;
t(1)   = 0;
z(1)   = -0.6;
y(1)   = -0.4;
h      = 0.1;
tfinal = 1;
m      = (tfinal-t(1))/h;
for j=1:m
    k11      = h*feval( f1 , t(j)     , y(j)     , z(j)    );
    k21      = h*feval( f2 , t(j)     , y(j)     , z(j)    );
    k12      = h*feval( f1 , t(j)+h/2 , y(j)+k11/2 , z(j)+k21/2);
    k22      = h*feval( f2 , t(j)+h/2 , y(j)+k11/2 , z(j)+k21/2);
    k13      = h*feval( f1 , t(j)+h/2 , y(j)+k12/2 , z(j)+k22/2);
    k23      = h*feval( f2 , t(j)+h/2 , y(j)+k12/2 , z(j)+k22/2);
    k14      = h*feval( f1 , t(j)+h   , y(j)+k13   , z(j)+k23  );
    k24      = h*feval( f2 , t(j)+h   , y(j)+k13   , z(j)+k23  );
    y(j+1)   = y(j) + (k11 + 2*k12 + 2*k13 + k14)/6;
    z(j+1)   = z(j) + (k21 + 2*k22 + 2*k23 + k24)/6;
    t(j+1)   = t(j) + h;
end
figure(1)
plot(t,[y;z])

fun = @(t,y)[y(2);sin(t)*exp(2*t) - 2*y(1) + 2*y(2)];
tspan = 0:h:tfinal;
[T,Y] = ode15s(fun,tspan,[y(1);z(1)],odeset('RelTol',1e-8,'AbsTol',1e-8));
figure(2)
plot(T,abs(Y-[y.',z.']))

[tahar zehrouri]

thank you very much bro

you saved me