Posting here vs math.stackexchange because I think my issue is syntax:
I'm trying to analyze the 2nd order ODE: y'' + 2y' + 2y = e^(-x) * sin(x) using MATLAB code for the midpoint method. I first converted the ODE to a system of 1st order equations and then tried to apply it below, but as the discretizations [m] are increased, the output is stopping at .2718. For example, m=11 yields:
ans =
0.2724
and m=101:
ans =
0.2718
and m=10001
ans =
0.2718
Here's the code:
function [y,t] = ODEsolver_midpointND(F,y0,a,b,m)
if nargin < 5, m = 11; end
if nargin < 4, a = 0; b = 1; end
if nargin < 3, a = 0; b = 1; end
if nargin < 2, error('invalid number of inputs'); end
t = linspace(a,b,m)';
h = t(2)-t(1);
n = length(y0);
y = zeros(m,n);
y(1,:) = y0;
for i=2:m
Fty = feval(F,t(i-1),y(i-1,:));
th = t(i-1)+h/2;
y(i,:) = y(i-1,:) + ...
h*feval(F,th,y(i-1,:)+(h/2)*Fty );
end
Separate file:
function F = Fexample1(t,y)
F1 = y(2);
F2 = exp(-t).*sin(t)-2.*y(2)-2.*y(1);
F = [F1,F2];
Third file:
[Y,t] = ODEsolver_midpointND('Fexample1',[0 0],0,1,11);
Ye = [(1./2).*exp(-t).*(sin(t)-t.*cos(t)) (1./2).*exp(-t).*((t-1).*sin(t)- t.*cos(t))];
norm(Y-Ye,inf)
Your ODE solver looks to me like it should work - however there's a typo in the analytic solution you're comparing to. It should be
Ye = [(1./2).*exp(-t).*(sin(t)-t.*cos(t)) (1./2).*exp(-t).*((t-1).*sin(t)+ t.*cos(t))];
i.e. with a +
sign before the t.*cos(t)
term in the derivative.