MATLAB TUTORIAL, part 2.4: Numerical Methods

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

  [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])