>> help ode45
 ode45  Solve non-stiff differential equations, medium order method.
    [TOUT,YOUT] = ode45(ODEFUN,TSPAN,Y0) with TSPAN = [T0 TFINAL] integrates 
    the system of differential equations y' = f(t,y) from time T0 to TFINAL 
    with initial conditions Y0. ODEFUN is a function handle. For a scalar T
    and a vector Y, ODEFUN(T,Y) must return a column vector corresponding 
    to f(t,y). Each row in the solution array YOUT corresponds to a time 
    returned in the column vector TOUT.  To obtain solutions at specific 
    times T0,T1,...,TFINAL (all increasing or all decreasing), use TSPAN = 
    [T0 T1 ... TFINAL].     
    
    [TOUT,YOUT] = ode45(ODEFUN,TSPAN,Y0,OPTIONS) solves as above with default
    integration properties replaced by values in OPTIONS, an argument created
    with the ODESET function. See ODESET for details. Commonly used options 
    are scalar relative error tolerance 'RelTol' (1e-3 by default) and vector
    of absolute error tolerances 'AbsTol' (all components 1e-6 by default).
    If certain components of the solution must be non-negative, use
    ODESET to set the 'NonNegative' property to the indices of these
    components.
    
    ode45 can solve problems M(t,y)*y' = f(t,y) with mass matrix M that is
    nonsingular. Use ODESET to set the 'Mass' property to a function handle 
    MASS if MASS(T,Y) returns the value of the mass matrix. If the mass matrix 
    is constant, the matrix can be used as the value of the 'Mass' option. If
    the mass matrix does not depend on the state variable Y and the function
    MASS is to be called with one input argument T, set 'MStateDependence' to
    'none'. ODE15S and ODE23T can solve problems with singular mass matrices.  
 
    [TOUT,YOUT,TE,YE,IE] = ode45(ODEFUN,TSPAN,Y0,OPTIONS) with the 'Events'
    property in OPTIONS set to a function handle EVENTS, solves as above 
    while also finding where functions of (T,Y), called event functions, 
    are zero. For each function you specify whether the integration is 
    to terminate at a zero and whether the direction of the zero crossing 
    matters. These are the three column vectors returned by EVENTS: 
    [VALUE,ISTERMINAL,DIRECTION] = EVENTS(T,Y). For the I-th event function: 
    VALUE(I) is the value of the function, ISTERMINAL(I)=1 if the integration 
    is to terminate at a zero of this event function and 0 otherwise. 
    DIRECTION(I)=0 if all zeros are to be computed (the default), +1 if only 
    zeros where the event function is increasing, and -1 if only zeros where 
    the event function is decreasing. Output TE is a column vector of times 
    at which events occur. Rows of YE are the corresponding solutions, and 
    indices in vector IE specify which event occurred.    
 
    SOL = ode45(ODEFUN,[T0 TFINAL],Y0...) returns a structure that can be
    used with DEVAL to evaluate the solution or its first derivative at 
    any point between T0 and TFINAL. The steps chosen by ode45 are returned 
    in a row vector SOL.x.  For each I, the column SOL.y(:,I) contains 
    the solution at SOL.x(I). If events were detected, SOL.xe is a row vector 
    of points at which events occurred. Columns of SOL.ye are the corresponding 
    solutions, and indices in vector SOL.ie specify which event occurred. 
 
    Example    
          [t,y]=ode45(@vdp1,[0 20],[2 0]);   
          plot(t,y(:,1));
      solves the system y' = vdp1(t,y), using the default relative error
      tolerance 1e-3 and the default absolute tolerance of 1e-6 for each
      component, and plots the first component of the solution. 
    
    Class support for inputs TSPAN, Y0, and the result of ODEFUN(T,Y):
      float: double, single
 
    See also ode23, ode113, ode15s, ode23s, ode23t, ode23tb, ode15i,
             odeset, odeplot, odephas2, odephas3, odeprint, deval,
             odeexamples, rigidode, ballode, orbitode, function_handle.

    Documentation for ode45

help vdp1
 vdp1  Evaluate the van der Pol ODEs for mu = 1
 
    See also ode113, ode23, ode45.

help vdp100
--- vdp100 not found. Showing help for vdp1000 instead. ---

 vdp1000  Evaluate the van der Pol ODEs for mu = 1000.
 
    See also ode15s, ode23s, ode23t, ode23tb.

help vdp1000
 vdp1000  Evaluate the van der Pol ODEs for mu = 1000.
 
    See also ode15s, ode23s, ode23t, ode23tb.

help vdpode
 vdpode  Parameterizable van der Pol equation (stiff for large MU).
    For the default value of MU = 1000 the equation is in relaxation
    oscillation, and the problem becomes very stiff.  The limit cycle has
    portions where the solution components change slowly and the problem is
    quite stiff, alternating with regions of very sharp change where it is
    not stiff (quasi-discontinuities). The initial conditions are close to
    an area of slow change so as to test schemes for the selection of the
    initial step size.
 
    The nested function J(T,Y) returns the Jacobian matrix dF/dY evaluated
    analytically at (T,Y). By default, the stiff solvers of the ODE Suite
    approximate Jacobian matrices numerically. However, if the ODE Solver
    property Jacobian is set to @J with ODESET, a solver calls the function
    to obtain dF/dY. Providing the solvers with an analytic Jacobian is not
    necessary, but it can improve the reliability and efficiency of
    integration.
 
    L. F. Shampine, Evaluation of a test set for stiff ODE solvers, ACM
    Trans. Math. Soft., 7 (1981) pp. 409-420.
 
    See also ode15s, ode23s, ode23t, ode23tb, odeset, function_handle.

[t,y]=ode45(@vdp1,[0 20],[2 0]);
plot(t,y(:,1));
Warning: MATLAB has disabled some advanced graphics rendering features by
switching to software OpenGL. For more information, click here. 
figure
plot(t,y(:,2));
[t,y]=ode45(@vdp1000,[0 20],[2 0]);
plot(t,y(:,1));
plot(t,y(:,2));
plot(t(1:00),y(1:100,1));
Error using plot
Vectors must be the same length. 
plot(t(1:100),y(1:100,1));
plot(t,y(:,1));
[t,y]=ode45(@vdp1,[0 6*pi],[2 0]);
plot(t,y(:,1));
[t,y]=ode45(@vdp1000,[0 4000*pi],[2 0]);
Operation terminated by user during ode45 (line 451) 

[t,y]=ode45(@vdp1000,[0 1000],[2 0]);
plot(t,y(:,1));
[t,y]=ode45(@vdp1000,[0 2000],[2 0]);
plot(t,y(:,1));
figure
plot(t,y(:,2));
plot(t(1:1000000),y(1:1000000,2));
plot(t(1:2000000),y(1:2000000,2));
plot(t(1:1500000),y(1:1500000,2));
[t,y]=ode23s(@vdp1000,[0 2000],[2 0]);
plot(t,y(:,1));
[t1,y1]=ode15s(@vdp1000,[0 2000],[2 0]);
plot(t1,y1(:,1));
plot(t1,y1(:,1),t,y(:,1));
[t1,y1]=ode15s(@vdp1000,t,[2 0]);
plot(t1,y1(:,1)-y(:,1));
[t,y]=ode23s(@vdp1000,[0 6000],[2 0]);
[t1,y1]=ode15s(@vdp1000,t,[2 0]);
plot(t1,y1(:,1)-y(:,1));
plot(t1,y1(:,1),t,y(:,1));
A=1;B=1;
F=@(t,Y)[A+Y(1).^2.*Y(2)-(B+1)*Y(1);B*Y(1)-Y(1).^2.*Y(2)]
F =
  function_handle with value:
    @(t,Y)[A+Y(1).^2.*Y(2)-(B+1)*Y(1);B*Y(1)-Y(1).^2.*Y(2)]
[t,y]=ode45(F,[0,20],[1;1]);
plot(t1,y1(:,1));
plot(t,y(:,1));
plot(t,y(:,2));
B=5;
F=@(t,Y)[A+Y(1).^2.*Y(2)-(B+1)*Y(1);B*Y(1)-Y(1).^2.*Y(2)]
F =
  function_handle with value:
    @(t,Y)[A+Y(1).^2.*Y(2)-(B+1)*Y(1);B*Y(1)-Y(1).^2.*Y(2)]
[t,y]=ode45(F,[0,20],[1;1]);
clear t1 y1
plot(t,y(:,1));
figure
plot(t,y(:,2));
figure
plot(y(:,1),y(:,2));
[t,y]=ode45(F,[0,40],[1;1]);
plot(y(:,1),y(:,2));
B=10;
F=@(t,Y)[A+Y(1).^2.*Y(2)-(B+1)*Y(1);B*Y(1)-Y(1).^2.*Y(2)]
F =
  function_handle with value:
    @(t,Y)[A+Y(1).^2.*Y(2)-(B+1)*Y(1);B*Y(1)-Y(1).^2.*Y(2)]
[t,y]=ode45(F,[0,40],[1;1]);
plot(y(:,1),y(:,2));
B=20;
F=@(t,Y)[A+Y(1).^2.*Y(2)-(B+1)*Y(1);B*Y(1)-Y(1).^2.*Y(2)]
F =
  function_handle with value:
    @(t,Y)[A+Y(1).^2.*Y(2)-(B+1)*Y(1);B*Y(1)-Y(1).^2.*Y(2)]
[t,y]=ode45(F,[0,40],[1;1]);
plot(y(:,1),y(:,2));
plot(t,y(:,1));
plot(t,y(:,2));
[t,y]=ode45(F,[0,200],[1;1]);
plot(y(:,1),y(:,2));
[t,y]=ode23s(F,[0,200],[1;1]);
plot(y(:,1),y(:,2));
B=1.6;
F=@(t,Y)[A+Y(1).^2.*Y(2)-(B+1)*Y(1);B*Y(1)-Y(1).^2.*Y(2)];
[t,y]=ode45(F,[0,200],[1;1]);
plot(y(:,1),y(:,2));
B=1.9;
F=@(t,Y)[A+Y(1).^2.*Y(2)-(B+1)*Y(1);B*Y(1)-Y(1).^2.*Y(2)];
[t,y]=ode45(F,[0,200],[1;1]);
plot(y(:,1),y(:,2));
B=2;
F=@(t,Y)[A+Y(1).^2.*Y(2)-(B+1)*Y(1);B*Y(1)-Y(1).^2.*Y(2)];
[t,y]=ode45(F,[0,200],[1;1]);
plot(y(:,1),y(:,2));
[t,y]=ode23s(F,[0,200],[1;1]);
plot(y(:,1),y(:,2));