 msSol  = tsolve(msSol , dT);
        refSol = tsolve(refSol, dT);

        <span class="comment">% Update the pressure based on the new saturation</span>
        <span class="comment">%msSol    = psolve(msSol);</span>
        <span class="comment">%refSol   = r_psolve(refSol);</span>

        <span class="comment">% store saturation at source</span>
        msSat = [msSat ; msSol.s(c2,1)];
        refSat = [refSat ; refSol.s(c2,1)];

        <span class="comment">% Increase time and continue if we do not w