% Mat_ODE.m

% This is a template file for solving a system of ordinary differential
% equations. This example is the van der Pol oscillator.

function [t,y] = call_RHS();
   tmax = 50;
   tspan = [0,tmax];
   y10=2; y20=1;  % initial conditions
   y0 = [y10 y20];   % initial condition vector
   [t,y] = ode23(@RHS,tspan,y0);

%  ode23 is an explicit Runge-Kutta (2,3) implementation for nonstiff ODEs
%  ode45 is a variable-stepsize implementation for nonstiff ODEs
%  ode15s is a variable-stepsize implementation for nonstiff ODEs
%  ode113 is a high-order variable-step Adams-Bashforth-Moulton solver
%     of orders 1 to 13.

%  change defaults for plotting
   set(0,'DefaultAxesFontSize',18)

%  set up variable vectors
   y1=y(:,1); y2=y(:,2);

%  Time course figure
   figure;

%  Use a 2-panel figure with 1 row and 2 columns
%  Left panel
   subplot(1,2,1);
   h=plot(t,y1,'black');
   set(h,'LineWidth',2);
   axis([0 tmax -3 3]);
   box off;
   xlabel('t');
   ylabel('y1');
  
%  Right panel
   subplot(1,2,2);
   h=plot(t,y2,'black');
   set(h,'LineWidth',2);
   axis([0 tmax -3 3]);
   box off;
   xlabel('t');
   ylabel('y2');

%  Define the RHS function for the ODEs
   function dydt = RHS(t,y);
 
%     set up the variable vectors
      y1=y(1); y2=y(2);

%     Parameters
      mu=1;

%     Differential equations
      y1dot = y2;
      y2dot = mu*(1-y1^2)*y2-y1;
      dydt=[y1dot;y2dot];
   end
end 
