% y' + 2y = sinx, y(0) = 0
% y = (e^(-2x) + 2sinx - cosx)/5

ya = @(x) (exp(-2*x) + 2*sin(x) - cos(x))/5;

vec = 0:0.001:1;
yaval = ya(vec);

yaval(1:10)

yi = zeros(1001,1001);
sumyi = zeros(1001,1);

[xi, yi] = ode45(@right3, [0:0.001:1], zeros(1001,1));

for i=1:1001
    for j=1:1001
        sumyi(j,1) = sumyi(j,1) + yi(i,j);
    end

    diff = yaval-sumyi;

    if max(diff(:)) < 0.01
        i
        break;
    else
        %max(diff(:))
    end
end

----

function [ dyi ] = right3(x, yi)

dyi = zeros(1001, 1);
for i = 0:1000
    dyi(i+1) = -2*yi(i+1) + (-1)^i * x^(2*i+1) / factorial(2*i+1);
end

end