function [Ux_All, Uy_All, M_All_b, M_All_u, Time_All, PosX, PosY, Front_All, FrontVelocity_All, Global_Area_All,Area_All,InitialMass,str] = GetSolutionStokesBFMyosin(config,SelectedTime)

global V0 M0 K A0 FX FY Actin_Velocity FM NormalX NormalY TotalArea  DomainX DomainY AreaExponent;

ticID = tic;

V0                   = config.V0;
M0                   = config.M0;
K                    = config.K;
A0                   = config.A0;
Domain               = config.Domain;
D                    = config.D; 
k_off                = config.k_off;
k_on                 = config.k_on;
eta                   = config.eta;                        
ksi                  = config.ksi;                                          
InitCondition        = config.InitCondition;
T_max                = config.T_max; 
MeshSize             = config.MeshSize;
Front.Formula        = config.FrontFormula;
VelBoundaryCond      = config.VelBoundaryCond;
NumFrames            = config.NumFrames;
ExactPolygon         = [];
sigma                   = config.sigma;
MaxMyosinAdvection   = config.MaxMyosinAdvection;
MaxMyosinDiffusion   = config.MaxMyosinDiffusion;
AreaExponent         = config.AreaExponent;
MyosinBoundaryCond   = config.MyosinBoundaryCond;

DomainX = Domain{1};
DomainY = Domain{2};
NumNodesX = round((DomainX(2)-DomainX(1))/MeshSize);
NumNodesY = round((DomainY(2)-DomainY(1))/MeshSize);
NumNodes  = NumNodesX*NumNodesY;
PosX = repmat((1:NumNodesX)', 1, NumNodesY);
PosX = DomainX(1)+(PosX(:) - 1/2)*MeshSize;
PosY = repmat(1:NumNodesY, NumNodesX,1);
PosY = DomainY(1)+(PosY(:) - 1/2)*MeshSize;

% Set time step size  
MaxVelocity = 1;
dt = config.TimeStepScale*MeshSize/(MaxVelocity*10);
NumTimeNodes = round(T_max/dt);
dt = T_max/NumTimeNodes; 
NumTimeNodes=round(T_max/dt);


% store the number of nonlinear solver iterations (segregated solver)
SegSolverIterations = zeros(NumTimeNodes,1);

% modify for precise last time output
det_last_output = 0;
iT_last_output = 0;
if  mod(NumTimeNodes,max(round(NumTimeNodes/NumFrames),1))>0
  iT_last_output = NumTimeNodes;
  det_last_output = 1;
  NumFramesFinal = floor(NumTimeNodes/max(round(NumTimeNodes/NumFrames),1))+2;
else
  NumFramesFinal = NumTimeNodes/max(round(NumTimeNodes/NumFrames),1)+1;
end

% file name
if config.FixedBoundary
    str = [config.VelBoundaryCond '_' config.MyosinBoundaryCond '_FxdFrt_N=' num2str(config.NumNodesX) '_T=' ...
        num2str(config.T_max) '_v0=' num2str(config.v0) '_mu_tot=' num2str(config.mu) ...
        '_alpha=' num2str(config.alpha) '_kon=' num2str(config.k_on) ...
        '_koff=' num2str(config.k_off)];    
else
    str = [config.VelBoundaryCond '_' config.MyosinBoundaryCond '_N=' num2str(config.NumNodesX) '_T=' ...
        num2str(config.T_max) '_v0=' num2str(config.v0) '_mu_tot=' num2str(config.mu) ...
        '_alpha=' num2str(config.alpha) '_kon=' num2str(config.k_on) ...
        '_koff=' num2str(config.k_off)];
end

if config.AvoidSingbyDiffu 
 str = [str '_dif-' num2str(MaxMyosinDiffusion)];
end

if config.AvoidSingbyAdvec
 str = [str '_adv-' num2str(MaxMyosinAdvection)];
end

if config.SolveMyosinsInOneStep
    str = [str '_cpld']; 
else
    str = [str '_segr'];
end

CurrentDir = pwd;
if ispc
    delim_ch = '\'; 
else
    delim_ch = '/'; 
end
if ~exist('myosin', 'dir')
    mkdir([CurrentDir delim_ch 'myosin' delim_ch]);
end

% preparing front velocity at t=0
if strcmp(config.VelBoundaryCond, 'ZeroDirichlet')
  FM = InitCondition.M;
else
  if config.bilinearInterp
    Actin_Velocity = @(x,y) [k_on/(k_on+k_off)*InitCondition.Ux(x,y) k_on/(k_on+k_off)*InitCondition.Uy(x,y)];          
  else
    FX = @(x,y) k_on/(k_on+k_off)*InitCondition.Ux(x,y);
    FY = @(x,y) k_on/(k_on+k_off)*InitCondition.Uy(x,y);
  end
end

switch config.InitDomainShape
  case 'Circle'
    TotalArea = config.InitArea;
    NormalX   = @(x,y) x/norm([x y]);
    NormalY   = @(x,y) y/norm([x y]);
  case 'ShapeFromSavedData'
    TotalArea = config.InitArea;
    iNeigh = 1;
    InitFront = config.InitFront;
    Tangents  = [InitFront.Polygon(end-iNeigh+1:end,:); InitFront.Polygon(1:end-iNeigh,:)] - [InitFront.Polygon(1+iNeigh:end,:); InitFront.Polygon(1:iNeigh,:)];
    Normals   = [-Tangents(:,2), Tangents(:,1)]./repmat(sqrt(Tangents(:,1).^2+Tangents(:,2).^2),1,2);
    NormalX   = scatteredInterpolant(InitFront.Polygon(:,1), InitFront.Polygon(:,2), Normals(:,1),'nearest');
    NormalY   = scatteredInterpolant(InitFront.Polygon(:,1), InitFront.Polygon(:,2), Normals(:,2),'nearest'); 
end


% Initialize front, mesh & solution     

[f, Front] = InitializeFront(Front, ExactPolygon, Domain, ...
  round(config.FrontRefinementScale*NumNodesX), T_max, config.FrontPropagate, ...
  config.FrontRedistMethod, config.FrontRedistFrequency);
if isempty(Front.Polygon)
    disp('Front variable is empty. Please choose a smaller MeshSize.')
    Ux_All=[]; Uy_All=[]; M_All_b=[]; M_All_u=[]; Time_All=[]; PosX=[]; PosY=[]; Front_All=[]; FrontVelocity_All=[]; Global_Area_All=[];
    return;
end

[Ux, Uy, M_b, M_u, InsideOldMask, InteriorOldMask, BoundaryOldMask] =...
    InitializeSolutionBFMyosin(PosX, PosY, Front, InitCondition, NumNodesX);
[NatNeighborPairs, NatNeighborBoundaryPairs, Vertices, VoronoiCells, Global_Area, D_Inv, Directions, MidPoints] =...
    InitializeControlVolumes(PosX, PosY, InsideOldMask, InteriorOldMask, BoundaryOldMask, Front, MeshSize, NumNodesX, NumNodesY);

InitialMass = sum(M_b.*Global_Area) + sum(M_u.*Global_Area);
% making sure that the initial solution has the correct M_total
if strcmp(config.InitialCondType,'SolutionFromSavedData')
  M_b = config.M_tot/InitialMass*M_b;
  M_u = config.M_tot/InitialMass*M_u;
end
  
% modify data structures for precise last time output
if config.returnAllData
    Ux_All = sparse(NumNodes, NumFramesFinal);
    Ux_All(:,1) = Ux;
    Uy_All = sparse(NumNodes, NumFramesFinal);
    Uy_All(:,1) = Uy;
    M_All_b = sparse(NumNodes, NumFramesFinal);
    M_All_b(:,1) = M_b;
    M_All_u = sparse(NumNodes, NumFramesFinal);
    M_All_u(:,1) = M_u;    
    Global_Area_All = sparse(NumNodes, NumFramesFinal);
    Global_Area_All(:,1) = Global_Area;
    Time_All = zeros(1, NumFramesFinal);
    Time_All(1,1) = SelectedTime;
    Area_All = zeros(1, NumFramesFinal);
    Area_All(1) = sum(Global_Area(InsideOldMask));
    Front_All = cell(1, NumFramesFinal);
    Front_All{1} = Front;
    FrontVelocity_All = cell(1, NumFramesFinal);
    FrontVelocity = zeros(size(Front.Polygon));
    for iVertex = 1:size(Front.Polygon,1)
        FrontVelocity(iVertex,:) = FT_velo_func(0, Front.Polygon(iVertex,:));
    end
    FrontVelocity_All{1} = FrontVelocity;
end

A_old = Global_Area(InsideOldMask);

fprintf('\n')
strlength = 0;
doVoronoi = false;
iNeigh = 1;

% segregated solver divergence detection
DivergenceCount = 0;

for iT = 1:NumTimeNodes
  
    if DivergenceCount > 10
      disp('Breaking the time iteration loop.');
      fprintf('Segregated solver diverged in the past %d iterations!',DivergenceCount);
       break;  
    end
    
    time = dt*iT;
    if ~mod(iT,round(NumTimeNodes/1000)) && ~config.DisplayFixedPointError
        strlength = DisplayProgress(ticID, iT, NumTimeNodes, strlength);
    end
    
    if ~config.FixedBoundary

        TotalArea = sum(A_old);
        if config.smoothing
            iNeigh = max(round(0.5*sqrt(size(Front.Polygon,1))),1);
        end

        Tangents = [Front.Polygon(end-iNeigh+1:end,:); Front.Polygon(1:end-iNeigh,:)] - [Front.Polygon(1+iNeigh:end,:); Front.Polygon(1:iNeigh,:)];
        Normals = [-Tangents(:,2), Tangents(:,1)]./repmat(sqrt(Tangents(:,1).^2+Tangents(:,2).^2),1,2);
        NormalX = scatteredInterpolant(Front.Polygon(:,1), Front.Polygon(:,2), Normals(:,1),'nearest');
        NormalY = scatteredInterpolant(Front.Polygon(:,1), Front.Polygon(:,2), Normals(:,2),'nearest');
        if config.bilinearInterp
            switch VelBoundaryCond
                case 'ZeroStress'
                    Actin_Velocity = @(x,y) bilinear2D_vec(PosX, PosY, k_on/(k_on+k_off)*Ux, k_on/(k_on+k_off)*Uy, InteriorOldMask, x, y, MeshSize);
                case 'ZeroDirichlet'
                    FM = @(x,y) bilinear2D(PosX, PosY, M_b+M_u, InteriorOldMask, x, y, MeshSize); 
            end
        else
            switch VelBoundaryCond
                case 'ZeroStress'
                    FX = scatteredInterpolant(PosX(InsideOldMask), PosY(InsideOldMask), k_on/(k_on+k_off)*Ux(InsideOldMask),'nearest');
                    FY = scatteredInterpolant(PosX(InsideOldMask), PosY(InsideOldMask), k_on/(k_on+k_off)*Uy(InsideOldMask),'nearest');
                case 'ZeroDirichlet'
                    FM = scatteredInterpolant(PosX(InsideOldMask), PosY(InsideOldMask), M_b(InsideOldMask)+M_u(InsideOldMask),'nearest');
            end
        end

    end
        
    
% Propogate front
    if ~config.FixedBoundary || iT==1
        
        [Front, FrontHitsOuterBox] = PropogateFront(f, time, Front, ExactPolygon, Domain, MeshSize);
        if FrontHitsOuterBox, break; end

    % Identify inside nodes    

        [InsideMask, InteriorMask, BoundaryMask] = IdentifyNodeTypes(PosX, PosY, Front, NumNodesX);

    % Update front movement metrics    
        [InteriorCandidatesMask, InteriorCandidatesIdx, BoundaryCandidatesIdx,new_nodes,gone_nodes] = ...
            FrontMovementUpdates(InsideMask, InteriorMask, BoundaryMask, InsideOldMask, InteriorOldMask, BoundaryOldMask);

    % Voronoi cell decomposition (Old nodes, new front)
    % Natural neighbors array: [i j], voronoi control volumes: A_i, edges: S_ij, inverse distances between nat.neighbors: D_Inv_ij 

        if doVoronoi || ~isempty(gone_nodes) || ~isempty(new_nodes)

            % update the discretization metrics that are dependent on voronoi
            % decomposition
            [NatNeighborBoundaryPairs, Vertices, VoronoiCells, D_Inv, Directions, MidPoints, NatNeighborPairs, Global_Area] = ...
                UpdateDiscretizationMetrics(BoundaryCandidatesIdx, PosX, PosY, InsideOldMask, NumNodesX, NumNodesY, InteriorCandidatesIdx, MeshSize, InteriorCandidatesMask);

        end
        % update the discretization metrics that are dependent on front
        % movement.
        [S, Global_Area, BoundaryIdxGlobal,BoundaryCurveLengths] = UpdateAreaEdgeLength(BoundaryCandidatesIdx, MeshSize, ...
            VoronoiCells, Vertices, NatNeighborBoundaryPairs, Directions, MidPoints, Front,true, true, InteriorCandidatesIdx, Global_Area);    


    % "Local" enumeration of the inside nodes 
        [I,J,NumInsideNodes] = LocalEnumeration(InsideOldMask,NumNodes,NatNeighborPairs);

    end
    
% Generate matrices for discretized system of equations       
% Update concentrations and velocities at the old nodes

    A = Global_Area(InsideOldMask);
    Ux_local = Ux(InsideOldMask);
    Uy_local = Uy(InsideOldMask);
    M_local_b = M_b(InsideOldMask);
    M_local_u = M_u(InsideOldMask);
    
    % the viscous term
    [Mat_Laplacian,~] = DiffusionTerm(eta, I, J, S, D_Inv, NumInsideNodes,true);
    CellVariablekisi =  CreateCellVariable(ksi,NumInsideNodes);
    % the cell-substrate adhesion term
    [Mat_source] = SourceTerm(A, CellVariablekisi); 
    AU = Mat_Laplacian + Mat_source;    
    
    % the active stress term
    [~,vector_x,~,vector_y] = GradientTerm(NatNeighborPairs, PosX, PosY, I, J, S, NumInsideNodes,M_local_b);
    
    % system matrix and right-hand-side vectors of the velocity equations
    BUx = sigma*vector_x;
    BUy = sigma*vector_y;
    
    % apply the zero velocity boundary condition
    if strcmp(VelBoundaryCond, 'ZeroDirichlet')
      [AU,BoundaryIdxLocal] = ApplyDirichletBndCondMat(AU,InsideOldMask,NumInsideNodes,size(Ux),BoundaryIdxGlobal);
      BUx = ApplyDirichletBndCondVec(BUx,0,BoundaryIdxLocal);
      BUy = ApplyDirichletBndCondVec(BUy,0,BoundaryIdxLocal);
    end
       
    % book-keeping for the old time step values
    M_local_b_PTS = M_local_b;
    M_local_u_PTS = M_local_u;
    
    % error printing initiation
    if config.DisplayFixedPointError
        if config.FixedPointStopCrite
            disp(['===================t=' num2str(SelectedTime+time) '==================='])
            fprintf('iter.\tUx\t\t\t     Uy\t\t\t    M\n');
        end
    end

    
   % the segregated solver iteration loop
   for ss=1:config.FixedPointIterationsMax
       
      if config.FixedPointStopCrite
          % store the old values for error calculation
          M_local_b_old = M_local_b;
          M_local_u_old = M_local_u;
          Ux_local_old = Ux_local;
          Uy_local_old = Uy_local; 
      end
      
      % update velocities
      Ux_local = AU\BUx;
      Uy_local = AU\BUy;
      
      % update advection velocities for bound Myosin
      if config.AvoidSingbyAdvec
        detM_local = M_local_b<MaxMyosinAdvection;
        AdvectionField = {@(ux,uy,m) ux.*detM_local.*(1-m/MaxMyosinAdvection) @(ux,uy,m) uy.*detM_local.*(1-m/MaxMyosinAdvection)};
        EdgeVelocities = AdvectionVelocity(AdvectionField, NatNeighborPairs, PosX, PosY, I, J, NumInsideNodes, 'CentralDifference',Ux_local,Uy_local,M_local_b);
      else
        AdvectionField = {@(ux,uy) ux @(ux,uy) uy};
        EdgeVelocities = AdvectionVelocity(AdvectionField, NatNeighborPairs, PosX, PosY, I, J, NumInsideNodes, 'CentralDifference',Ux_local,Uy_local);
      end
      [Mat_advection,~] = AdvectionTerm(EdgeVelocities, I, J, S, NumInsideNodes); 
      
      
      % system matrix for the diffusion of unbound myosin 
      if config.AvoidSingbyDiffu   % non-linear diffusion coefficient
        detM_local = M_local_u<MaxMyosinDiffusion;
        DiffCoeff = @(m) D*detM_local.*(1-m/MaxMyosinDiffusion);
        [D_face] = DiffusionCoefficient(DiffCoeff,I,J,'HarmonicAverage',M_local_u);
        [Mat_diffusion,~] = DiffusionTerm(D_face, I, J, S, D_Inv, NumInsideNodes,false);

      else                         % constant diffusion coefficient
        
        [Mat_diffusion,~] = DiffusionTerm(D, I, J, S, D_Inv, NumInsideNodes,true);
          
      end
      
      % transient terms
      if ss == 1
        [Mat_transient,BM_b] = TransientTerm(A, A_old, M_local_b_PTS, true, NumInsideNodes);
        [~,BM_u] = TransientTerm(A, A_old, M_local_u_PTS, true, NumInsideNodes);
      end
      
      % reaction term (matrix, vector), myosin, and the solution      
      if config.SolveMyosinsInOneStep
          
          % the reaction matrices are the same as the transient term
          % matrix, so no need to build them again! 
          AM = [Mat_transient*(1+k_off*dt) + dt*(Mat_advection) , -k_on*dt*Mat_transient; ...
              -k_off*dt*Mat_transient , Mat_transient*(1+k_on*dt) + dt*(Mat_diffusion)];
          
          if strcmp(MyosinBoundaryCond,'ZeroFlux')
              M_local = AM\[BM_b;BM_u];
          else
              % correct for known fluxes from the membrane
              % the total flux is known: (M_b * u) * n | mem. = (M_b * u)*n
              BoundaryIdxLocal = GetBoundaryIndexLocal(InsideOldMask,size(M_local_b),BoundaryIdxGlobal);
              Flux = M_local_b(BoundaryIdxLocal).*...
                  (  NormalX(PosX(BoundaryIdxGlobal),PosY(BoundaryIdxGlobal)).*Ux_local(BoundaryIdxLocal) + ...
                  NormalY(PosX(BoundaryIdxGlobal),PosY(BoundaryIdxGlobal)).*Uy_local(BoundaryIdxLocal)  ).* ...
                  BoundaryCurveLengths;
              % add the flux to the RHS of advection equation
              RHS_b = BM_b;
              RHS_b(BoundaryIdxLocal) = RHS_b(BoundaryIdxLocal)-dt*Flux;
              
              % add the same flux to the RHS of diffusion equation
              RHS_u = BM_u;
              RHS_u(BoundaryIdxLocal) = RHS_u(BoundaryIdxLocal)+dt*Flux;
              
              % solve the system
              M_local = AM\[RHS_b;RHS_u];
          end
          
          M_local_b = M_local(1:NumInsideNodes);
          M_local_u = M_local(1+NumInsideNodes:end);
      else
          
          % solve bound and unbound myosin equations separately
          % the reaction are treated explicitly in both approaches
          % M_u ---k_on---> M_b
          [reactant1, product1] = ReactionVector({1,2,k_on}, A_old, M_local_u);
          % M_b ---k_off--> M_u
          [reactant2, product2] = ReactionVector({1,2,k_off}, A_old, M_local_b);            
          
          if strcmp(MyosinBoundaryCond,'ZeroFlux')          
          
              AM =  Mat_transient + dt*(Mat_diffusion);    
              M_local_u = AM\(BM_u + dt*(product2+reactant1));          

              AM = Mat_transient + dt*(Mat_advection);
              M_local_b = AM\(BM_b + dt*(product1+reactant2));
          
          else
              % correct for known fluxes from the membrane
              % the total flux is known: (M_b * u) * n | mem. = (M_b * u)*n
              BoundaryIdxLocal = GetBoundaryIndexLocal(InsideOldMask,size(M_local_b),BoundaryIdxGlobal);
              Flux = M_local_b(BoundaryIdxLocal).*...
                  (  NormalX(PosX(BoundaryIdxGlobal),PosY(BoundaryIdxGlobal)).*Ux_local(BoundaryIdxLocal) + ...
                  NormalY(PosX(BoundaryIdxGlobal),PosY(BoundaryIdxGlobal)).*Uy_local(BoundaryIdxLocal)  ).* ...
                  BoundaryCurveLengths;
              
              % diffusiona and transient term matrices
              AM =  Mat_transient + dt*(Mat_diffusion);
                        
              % add the flux to the RHS of diffusion equation
              RHS_u = BM_u + dt*(product2+reactant1);
              RHS_u(BoundaryIdxLocal) = RHS_u(BoundaryIdxLocal)+dt*Flux;              
              M_local_u = AM\RHS_u;          

              % advection and transient term matrices
              AM = Mat_transient + dt*(Mat_advection);
             
              % add the same flux to the RHS of advection equation
              RHS_b = BM_b + dt*(product1+reactant2);
              RHS_b(BoundaryIdxLocal) = RHS_b(BoundaryIdxLocal)-dt*Flux;
              M_local_b = AM\RHS_b;              
          end
          
      end
      
      % update the Myosin contractile stress
      [~,vector_x,~,vector_y] = GradientTerm(NatNeighborPairs, PosX, PosY, I, J, S, NumInsideNodes,M_local_b);
    
      
      % update the right-hand-side velocity vectors
      BUx = sigma*vector_x;
      BUy = sigma*vector_y;
      
      % apply the zero velocity boundary condition
      if strcmp(VelBoundaryCond, 'ZeroDirichlet')
        BUx = ApplyDirichletBndCondVec(BUx,0,BoundaryIdxLocal);
        BUy = ApplyDirichletBndCondVec(BUy,0,BoundaryIdxLocal);
      end
      
      % Error control and stopping criteria
      if config.FixedPointStopCrite                % Maximum absolute error
          
          % calculate solution errors
          MaxErrorUx = full(max(abs(Ux_local-Ux_local_old)));
          MaxErrorUy = full(max(abs(Uy_local-Uy_local_old)));
          MaxErrorMb = full(max(abs(M_local_b-M_local_b_old)));
          MaxErrorMu = full(max(abs(M_local_u-M_local_u_old)));
          MaxErrorM  = max(MaxErrorMb, MaxErrorMu);
          MaxErrorSol = max(max(MaxErrorUx,MaxErrorUy),MaxErrorM);
          if config.DisplayFixedPointError
              fprintf('%d\t%8.7e\t%8.7e\t%8.7e\n',ss,MaxErrorUx,MaxErrorUy,MaxErrorM); 
          end          
          if MaxErrorSol<config.FixedPointIterationsTol
              SegSolverIterations(iT) = ss;
             break;
          end
          if ss==config.FixedPointIterationsMax && MaxErrorSol>=config.FixedPointIterationsTol
              SegSolverIterations(iT) = ss;
              fprintf('Warning: the segregated solver did not converge at t = %5.3f! Max. Sol. error is : %8.7e after %i iterations\n',SelectedTime+time,MaxErrorSol,config.FixedPointIterationsMax);
              DivergenceCount = DivergenceCount+1;
          end
      end

            
   end % end segregated solver
    
     
    Ux(InsideOldMask) = Ux_local;
    Uy(InsideOldMask) = Uy_local;
    M_b(InsideOldMask) = M_local_b;
    M_u(InsideOldMask) = M_local_u;
    
% Adjust concentaration at the nodes neighboring those that disappear and new nodes
   if ~config.FixedBoundary
        [Ux, Uy, M_b, M_u, Global_Area, InsideOldMask, BoundaryOldMask, NatNeighborPairs] = GoneNodesAdjustmentBFMyosin(gone_nodes, Global_Area, Ux, Uy, M_b, M_u, PosX, PosY, InsideOldMask, BoundaryOldMask, NatNeighborPairs, NumNodesX, NumNodesY, MeshSize, Front);
        [Ux, Uy, M_b, M_u, Global_Area, InsideOldMask, NatNeighborPairs] = NewNodesAdjustmentBFMyosin(new_nodes, Global_Area, Ux, Uy, M_b, M_u, PosX, PosY, InsideOldMask, BoundaryOldMask, NatNeighborPairs, NumNodesX, NumNodesY, MeshSize, Front);
   end
   
% Collect the results

    if config.returnAllData && ~mod(iT,max(round(NumTimeNodes/NumFrames),1))
        ii = iT/max(round(NumTimeNodes/NumFrames),1)+1;
        Time_All(ii) = SelectedTime+time;
        Area_All(ii) = TotalArea;
        Front_All{ii} = Front;
        FrontVelocity = zeros(size(Front.Polygon));
        for iVertex = 1:size(Front.Polygon,1)
            FrontVelocity(iVertex,:) = FT_velo_func(time, Front.Polygon(iVertex,:));
        end
        FrontVelocity_All{ii} = FrontVelocity;
        Ux_All(:,ii) = Ux;
        Uy_All(:,ii) = Uy;
        M_All_b(:,ii) = M_b;
        M_All_u(:,ii) = M_u;
        Global_Area_All(:,ii) = Global_Area;
        
        if config.SaveDuringSolve
          
          % Save temporary
          dataMask = Time_All>0;
          dataMask(1) = true;

          M_All_temp_b = M_All_b;
          M_All_temp_u = M_All_u;
          Ux_All_temp = Ux_All;
          Uy_All_temp = Uy_All;
          Time_All_temp = Time_All;
          Area_All_temp = Area_All;
          Front_All_temp = Front_All;
          FrontVelocity_All_temp = FrontVelocity_All;
          Global_Area_All_temp = Global_Area_All;

          M_All_b = M_All_b(:,dataMask);
          M_All_u = M_All_u(:,dataMask);
          Ux_All = Ux_All(:,dataMask);
          Uy_All = Uy_All(:,dataMask);
          Time_All = Time_All(dataMask);
          Area_All = Area_All(dataMask);
          Front_All = Front_All(dataMask);
          FrontVelocity_All = FrontVelocity_All(dataMask);
          Global_Area_All = Global_Area_All(:,dataMask);      

          save([CurrentDir delim_ch 'myosin' delim_ch str '.mat'], 'config', 'M_All_b', 'M_All_u', 'Ux_All', 'Uy_All', 'Time_All', 'Area_All', ...
          'PosX', 'PosY', 'Front_All', 'FrontVelocity_All',...
          'Global_Area_All', 'config','InitialMass','SegSolverIterations','-v7.3', 'str');

          M_All_b = M_All_temp_b;
          M_All_u = M_All_temp_u;
          Ux_All = Ux_All_temp;
          Uy_All = Uy_All_temp;
          Time_All = Time_All_temp;
          Area_All = Area_All_temp;
          Front_All = Front_All_temp;
          FrontVelocity_All = FrontVelocity_All_temp;
          Global_Area_All = Global_Area_All_temp;

          clear M_All_temp_b M_All_temp_u Ux_All_temp Uy_All_temp Time_All_temp Area_All_temp Front_All_temp Global_Area_All_temp;
          
        end        
        
    end

    
    % modify for precise last time output    
    if config.returnAllData && (iT==iT_last_output) && det_last_output
        Time_All(ii+1) = SelectedTime+time;
        Area_All(ii+1) = TotalArea;
        Front_All{ii+1} = Front;
        FrontVelocity = zeros(size(Front.Polygon));
        for iVertex = 1:size(Front.Polygon,1)
            FrontVelocity(iVertex,:) = FT_velo_func(time, Front.Polygon(iVertex,:));
        end
        FrontVelocity_All{ii+1} = FrontVelocity;
        Ux_All(:,ii+1) = Ux;
        Uy_All(:,ii+1) = Uy;
        M_All_b(:,ii+1) = M_b;
        M_All_u(:,ii+1) = M_u;
        Global_Area_All(:,ii+1) = Global_Area;      
    end
    
% Visualize solution while running 
    
    if config.Visualization && ~mod(iT,round(NumTimeNodes/100))
        if iT == round(NumTimeNodes/100)
            figure('Position', [100 400 1400 500]);
        end
        subplot(1,2,1)
        VisualizeSolution(M_b, PosX, PosY, MeshSize, SelectedTime+time, Front, 0, 1)
        xlim(Domain{1})
        ylim(Domain{2})
        colorbar
        subplot(1,2,2)
        VisualizeSolution(sqrt(Ux.^2+Uy.^2), PosX, PosY, MeshSize,SelectedTime+time, Front, 0, 0)
        colorbar
        hold on
        VisualizeSolution({Ux, Uy, 1}, PosX, PosY, MeshSize, SelectedTime+time, Front, 1,0,{'color','w','LineWidth',1})
        xlim(Domain{1})
        ylim(Domain{2})    
        hold off
        pause(0.1)
    end
    
    BoundaryOldMask = BoundaryMask;
    InteriorOldMask = InteriorMask;
    A_old = Global_Area(InsideMask);
    if isempty(gone_nodes) && isempty(new_nodes)
        doVoronoi = false;
    else
        doVoronoi = true;
    end
end;
fprintf('\n')

if ~config.returnAllData
    Ux_All = Ux;
    Uy_All = Uy;
    M_All_b = M_b;
    M_All_u = M_u;
    Time_All = SelectedTime+time;
    Area_All = TotalArea;
    Front_All{1} = Front;
    Global_Area_All = Global_Area;
    FrontVelocity_All = zeros(size(Front.Polygon));
    for iVertex = 1:size(Front.Polygon,1)
        FrontVelocity_All(iVertex,:) = FT_velo_func(time, Front.Polygon(iVertex,:));
    end
end

if isempty(ExactPolygon)
  % free internal front representation
  fcmd.command = 'free';
  fcmd.front = f;
  fronTierCmd(fcmd);
end

dataMask = Time_All>0;
dataMask(1) = true;
M_All_b = M_All_b(:,dataMask);
M_All_u = M_All_u(:,dataMask);
Ux_All = Ux_All(:,dataMask);
Uy_All = Uy_All(:,dataMask);
Time_All = Time_All(dataMask);
Area_All = Area_All(dataMask);
Front_All = Front_All(dataMask);
Global_Area_All = Global_Area_All(:,dataMask);

% save the output file
save([CurrentDir delim_ch 'myosin' delim_ch str '.mat'], 'config', 'M_All_b', 'M_All_u', 'Ux_All', 'Uy_All', 'Time_All', 'Area_All', 'PosX', 'PosY', 'Front_All', 'FrontVelocity_All', 'Global_Area_All','InitialMass', 'config','str', 'SegSolverIterations','-v7.3');
str = [CurrentDir delim_ch 'myosin' delim_ch str ];
