Solution of (stiff) ordinary differential equations.
(t,y)=stiffeuler(func, t0, y0, h, nsteps)
(y)=stiffeuler(func, t0, y0, h, nsteps, tol)
func The function of the form (dybydt)=func(t, y) , that returns the derivative of the desired function.
t0 The initial time.
y0 Initial condition.
nsteps The number of steps.
tol The convergence criterion used for solution for nonlinear equations at each step.
t The vector of time instances at which the solution was found.
y The solution. Each row of y is the solution at the time instance in the corresponding row of t .

>>(t,y)=stiffeuler(@(t,y)[ y[2]; 10*(1-y[2]^2)*y[2]-y[1]] , 0,
>                    [1 0]',h=.25,round(50/h),1.e-5)
>>plot(t, y)
>>ylabel('y[1], y[2]')
>>title({'Implicit Euler Method for Stiff ODEs',"Van der Pol Oscillations"})