MATLAB EC
function [t,y] = euler(f,tspan,y0,N)
% Solves the IVP y' = f(t,y), y(t0) = y0 in the time interval tspan = [t0,tf]
% using Euler's method with N time steps.
% Input:
% f = name of inline function or function M-file that evaluates the ODE
% (if not an inline function, use: euler(@f,tspan,y0,N))
% For a system, the f must be given as column vector.
% tspan = [t0, tf] where t0 = initial time value and tf = final time value
% y0 = initial value of the dependent variable. If solving a system,
% initial conditions must be given as a vector.
% N = number of steps used.
% Output:
% t = vector of time values where the solution was computed
% y = vector of computed solution values.
m = length(y0);
t0 = tspan(1);
tf = tspan(2);
h = (tf-t0)/N; % evaluate the time step size
t = linspace(t0,tf,N+1); % create the vector of t values
y = zeros(m,N+1); % allocate memory for the output y
y(:,1) = y0'; % set initial condition
for n=1:N
y(:,n+1) = y(:,n) + h*f(t(n),y(:,n)); % implement Euler's method
end
t = t'; y = y'; % change t and y from row to column vectors
end
function [t,y] = impeuler(f,tspan,y0,N)
% Solves the IVP y? = f(t,y), y(t0) = y0 in the time interval tspan=[t0,tf]
% using improved Euler's method with N time steps.
% Input:
% f = name of inline function or function M-file that evaluates the ODE
% (if not an inline function, use: euler(@f,tspan,y0,N))
% For a system, the f must be given as column vector.
% tspan = [t0, tf] where t0 = initial time value and tf = final time value
% y0 = initial value of the dependent variable. If solving a system,
% initial conditions must be given as a vector.
% N = number of steps used. % Output:
% t = vector of time values where the solution was computed
% y = vector of computed solution values.
m = length(y0);
t0 = tspan(1);
tf = tspan(2);
h = (tf-t0)/N; % evaluate the time step size
t = linspace(t0,tf,N+1); % create the vector of t values
y = zeros(m,N+1); % allocate memory for the output y
y(:,1) = y0'; % set initial condition
for n=1:N
f1=f(t(n),y(:,n));
f2=f(t(n+1),y(:,n)+h*f1);
y(:,n+1) = y(:,n) + h/2*(f1+f2); % implement improved Euler's method
end
t = t'; y = y'; % change t and y from row to column vectors
end