title: Use intermediate value on ode MATLAB solver
body: <p>I am doing a time integration of a system of ODEs using a stiff solver (ode15s). It is working, but I want to speed things up.</p>
  
  <p>The system of equations is given in state space form:</p>
  
  <pre><code>function [dx] = fun(t,x,M,C,K,other_parameters)
      % Mx'' + Cx' + Kx = F(t)
      % BUNCH OF CALCULATIONS
      F = solveP(x,t);
  
      A = [zeros(n) eye(n) ; -M\K -M\C];
      b = M\F;
  
      dx = A*x + b
  end
  </code></pre>
  
  <p>The trick part here is the forcing function F. It is highly nonlinear and depends on the x and t parameters. It uses the x parameters to solve a Poisson-type 2D equation (by the Finite Volume method). The force F is proportional to the Poisson equation solution.</p>
  
  <pre><code>function [F] = solveP(x,t)
      % initialize solution
      Phi = zeros(Ni,Nj);
  
      % solve iteratively
      % ...
      while (~converged)
           % some calculations
  
           % iterative solver
           Phi(i,j) = (aE*Phi(i,j+1) + aW*Phi(i,j-1) + aN*Phi(i+1,j) +...
                       aS*Phi(i-1,j) + S(i,j))/aP;
      end
  
  
      % calculate F
      F = sum(Phi(:)); % discrete integration over domain
  end
  </code></pre>
  
  <p>Solving the Poisson equation by a iterative method requires an initial condition, which I set to zero (<code>Phi=zeros(Ni,Nj)</code>). I thought I could improve the speed of calculations by providing a better initial estimative of the ϕ field (a better initial condition would take faster to reach the sought answer). The optimal initial condition I can think (besides ϕ=0) is the value of the ϕ field obtained in the previous iteration (the last step) of the ode solver (<code>ϕ(k)_initial=ϕ(k-1)</code>).</p>
  
  <p>Bottom line is: how do I use/save intermediate values in the ode solution?</p>
  
  <p>PS: I tried using the persistent variables, but that is not a good solution. The ode solver calculates the function in several points before advancing in time. The persistent variable saves the converged ϕ field every time the ode calls the odefun fun. That is not exactly what I want and this actually provides wrong answers as time marching goes on.</p>