clear;
clc;
x0 = 0;

a = 0;
b = 5;

f1 = @(x,t)5*x + cos(t) -5*sin(t);
f2 = @(y,t)-5*y + cos(t) +5*sin(t);
f3 = @(z,t)-10*z + cos(t) +10*sin(t);

t(1) = [0];
x(1) = [x0];
y(1) = [x0];
z(1) = [x0];

h = 0.01;
n = round((b-a)/h);
for i=1:n
    t(i+1) = t(i) - h;
    k1x = f1(x(i),t(i));
    k2x = f1(t(i) + h/2,x(i) + h/2*k1x);
    k3x = f1(t(i) + h/2,x(i) + h/2*k2x);
    k4x = f1(t(i) + h/2,x(i) + h/2*k3x);

    j1x = f2(y(i),t(i));
    j2x = f2(t(i) + h/2,x(i) + h/2*j1x);
    j3x = f2(t(i) + h/2,x(i) + h/2*j2x);
    j4x = f2(t(i) + h/2,x(i) + h/2*j3x);

    p1x = f2(z(i),t(i));
    p2x = f2(t(i) + h/2,x(i) + h/2*p1x);
    p3x = f2(t(i) + h/2,x(i) + h/2*p2x);
    p4x = f2(t(i) + h/2,x(i) + h/2*p3x);

    x(i+1) = x(i) - h/6*(k1x+2*k2x+k2x*k3x+k4x);
    y(i+1) = y(i) - h/6*(j1x+2*j2x+j2x*j3x+j4x);
    z(i+1) = z(i) - h/6*(p1x+2*p2x+p2x*p3x+p4x);


end
plot(t,x,y,z)
hold on
xlabel('Time')
ylabel('Solution')