close all clear all T = 1; N=10; h = T/N; t = linspace(0,1,N+1)'; f = @(t,y)-y.^2+2*t.*y; y = @(t) exp(t.^2)./(1+sqrt(pi)*erfi(t)/2); %exact solution yFE = zeros(size(t));yFE(1)=1; %forward Euler yBE = zeros(size(t));yBE(1)=1; %backward Euler tFine = linspace(0,1,100)'; % fine mesh for plotting for n=1:N yFE(n+1) = yFE(n)+h*f(t(n), yFE(n)); yBE(n+1) = (2*t(n+1)-1/h + sqrt((2*t(n+1)-1/h)^2+4*yBE(n)/h))/2; % see book for reasoning plot(tFine,y(tFine),'b',t(1:n+1),yFE(1:n+1),'r-o', t(1:n+1), yBE(1:n+1),'k-x') axis([0 1 0.75 1.18]) xlabel('t','fontsize',14) ylabel('y','fontsize',14) legend({'Exakt losning', 'Euler framat', 'Euler bakat'},'location','NorthWest','fontsize',14) hold on pause() end %% Add direction field [tMesh, yMesh] = meshgrid(0:0.1:1., 0.75:0.02:1.2); F_t = ones(size(tMesh)); F_y = f(tMesh, yMesh); figure(2) quiver(tMesh, yMesh, F_t, F_y,1.,'b') axis tight; xlabel('t'), ylabel('y') title('Direction field for dy/dt = -y^2+2ty. That is (1,f(t,y))') hold on for n=1:N plot(tFine,y(tFine),'b',t(1:n+1),yFE(1:n+1),'r-o', t(1:n+1), yBE(1:n+1),'k-x','linewidth',2) axis([0 1 0.75 1.2]) xlabel('t','fontsize',14) ylabel('y','fontsize',14) legend({'(1,f(t,y))', 'Exakt losning', 'Euler framat', 'Euler bakat'},'location','NorthWest','fontsize',14) hold on pause() end