Runge-Kutta 2

So I'm doing a project for school, where i have to solve the x"+5x'+4x = 3 - 2t - t^2 equation with Runge-Kutta 2.

intervalmin = 0;
intervalmax = 1;
h = 0.1;
g = 0.02;
numnodes = ceil(((intervalmax - intervalmin) / h) + 1);
numnodes2 = ceil(((intervalmax - intervalmin) / g) + 1);
inival = 0;

t = zeros(1, numnodes);
x = zeros(2, numnodes);
t(1) = intervalmin;
x(:,1) = inival;
f = @(t,x) [x(2); -5 * x(2) - 4 * x(1) + 3 - 2 * (t) - t.^2];

for i = 2:numnodes
    t(i) = t(i - 1) + h;
    k1 = f(t(i - 1), x(:, i - 1));               
    k2 = f(t(i - 1) + h / 2, x(:, i-1) + (h / 2) *k1);
    k3 = f(t(i - 1) + h, x(:, i - 1) - h * k1 + 2 * h * k2);
    x(:, i) = x(:, i - 1) + (h / 6) * (k1 + 4 * k2 + k3); 
end

o = zeros(1, numnodes2);
z = zeros(2, numnodes2);
o(1) = intervalmin;
z(:, 1) = inival;
m = @(o, z) [z(2); -5 * z(2) - 4 * z(1) + 3 - 2 * (o) - o.^2];

for i = 2:numnodes2
    o(i) = o(i - 1) + g;
    k1 = m(o(i - 1), z(:, i - 1));               
    k2 = m(o(i - 1) + h / 2, z(:, i - 1) + (g / 2) * k1);
    k3 = m(o(i - 1) + h, z(:, i - 1) - g * k1 + 2 * g * k2);
    z(:, i)=z(:, i - 1) + (g / 6) * (k1 + 4 * k2 + k3);  
end

[TT, XX] = ode45(f,[intervalmin intervalmax],x(:, 1));
figure 
plot(t,x(2,:),'.-',o,z(2,:),'.-',TT,XX(:,2),'.-')
hold on
syms x(t)
dz = diff(x);
ode = diff(x,t,2) == (-5 * x(2) - 4 * x(1) + 3 - 2 * (t) - t.^2);
cond1 = x(0) == 1;
cond2 = dz(0) == 1;
dsolve(ode, cond1, cond2)

legend('h=0.1','h=0.02','ode45')
ax = gca;
ax.XAxisLocation = 'origin';
ax.YAxisLocation = 'origin';
grid on;

I'm pretty new to matlab, and just can't seem to find a way to make it work, so if somebody could help me solve it, i would be grateful

[Nagy Csaba Norbert]

One more question, i have a task with initialize conditions at x(0) = 1 and x'(0) = 1. Can someone also help with this?

[James Tursa]

This looks like a 3rd order method your are implementing. Is that what you are supposed to be doing? Or are you supposed to be implementing a 2nd order method as your title suggests?

Answer 1

[James Tursa]

Normally one would plot the position, not the velocity. So I would have expected you to plot the "1" index of your solutions. E.g.,

plot(t,x(1,:),'.-',o,z(1,:),'.-',TT,XX(:,1),'.-')