# MATLAB TUTORIAL, part 2.4: Numerical Methods

Solving ODEs

#### I. Verification

function solveDifferentialEquation
%This function solves y'+y = x from starttime to endtime with y(0) = intialCondition
%Parameter initializations
starttime = 0;
endtime = 20;
initialCondition = -2; function yprime = DifferentialEquationFunction(x,y)
yprime = y - x;
end[xVals,yVals] = ode45(@DifferentialEquationFunction,[starttime,endtime],initialCondition);
plot(xVals,yVals);
end %Function Complete


function solveDifferentialEquation2
% Function to solve the differential equation
% of a particle in motion due to impact of an external sinusoidal
% force ( yprime = F*sin(w*t) - C*y )
% from a starttime to stoptime, where the intial y value is assigned as
% initialCondition
%Parameter initializations
n = 5;
onesv = [1, 1.45, 1.34, 1.81, 3.5]; %acts as a scalar vector that scales the parameters by its %entries.
C = 10.*onesv;
F = 1.*onesv;
w = 1.*onesv; %With scaled frequencies,w, the plots should show a trend of increasing spread with t5he 3rd plot as an exception since 1.34 defies the pattern.

starttime = 0;
endtime = 20;
initialCondition = 0;
counter = 1;
%Loop to solve the ODE and plot five different graphs.
for indeX = 1:n
[t,y] = ode45(@DifferentialEquationFunction2,[starttime,endtime],initialCondition);
counter = counter + 1;
figure(indeX)
plot(t,y);
title('Solution to differential equation')
xlabel('time')
ylabel('y-value')
end
function yprime = DifferentialEquationFunction2(t,y)      % Function defining the ODE
yprime = F(counter)*sin(w(counter)*t) - C(counter)*y;
end
end          %Function Complete

Direction fields

[tm,ym]=meshgrid(tval,yval);
dt = tval(2) - tval(1);
dy = yval(2) - yval(1);
fv = vectorize(f);
if isa(f,'function_handle')
fv = eval(fv);
end
yp=feval(fv,tm,ym);
s = 1./max(1/dt,abs(yp)./dy)*0.35;
h = ishold;
quiver(tval,yval,s,s.*yp,0,'.r'); hold on;
quiver(tval,yval,-s,-s.*yp,0,'.r');
if h
hold on
else
hold off
end
axis([tval(1)-dt/2,tval(end)+dt/2,yval(1)-dy/2,yval(end)+dy/2])