set the reproducibility report with all the notebooks. The supplementary...

set the reproducibility report with all the notebooks. The supplementary material file must be updated

LICENSE

+
+MIT License
+
+Copyright (c) 2021 Davide Torlo, Philipp Öffner and Hendrik Ranocha
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.

README.md

 # Stability of Positivity Preserving Patankar-Type Schemes
 
-This repository contains the code for the reproducibility of "A new Stability Approach for Positivity-Preserving Patankar-type Schemes" by D. Torlo, P. Öffner and H. Ranocha.
+[![License: MIT](https://img.shields.io/badge/License-MIT-success.svg)](https://opensource.org/licenses/MIT)
\ No newline at end of file
+
+This repository contains the code for the reproducibility of "A new Stability Approach for Positivity-Preserving Patankar-type Schemes" by D. Torlo, P. Öffner and H. Ranocha [arXiv preprint](https://arxiv.org).
+
+We provide some *Julia* codes and Jupyter notebooks to reproduce the simulations of the paper, some *Mathematica* notebook to reproduce the symbolic computations and the results of the paper, and the supplementary material file with all the results that could not fit in the article.
+
+### Julia codes
+In folder [simulations](simulations/) there are:
+* the implementations of all the Patankar and modified Patankar schemes ([Patankar Runge--Kutta](simulations/mPRK.jl) and [modified Patankar Deferred Correction](simulations/mPDeC.jl)) and some auxiliary functions ([Lagrange polynomials file](simulations/Lagrange_polynomials.jl))
+* the study of the general **linear $2\times 2$ system** searching for the maximum $\Delta t$ bound for which we do not observe oscillations for all the schemes is in the [CFL notebook](simulations/SearchingCFLconditions.ipynb)
+* the study of the spurious steady state for all the methods is in the [Inconsistency notebook](simulations/Inconsistency.ipynb)
+* the test of the stability bound previously found on **scalar nonlinear problem** is in [Nonlinear Problem notebook](simulations/Nonlinear_Problems.ipynb)
+* the test of the inconsistency on **Robertson problem** is in [Robertson notebook](simulations/Robertson_Problem.ipynb)
+* the test of the inconsistency on **HIRES problem** is in [HIRES notebook](simulations/Problem_HIRES.ipynb)
+
+
+### Mathematica notebooks
+In folder [symbolic_computations](simbolic_computations/) contains:
+* the computations of the $\Delta t$ bound on the general **linear $2\times 2$ system** for the **MPRK(2,2,1)** scheme is in the [MPRK(2,2,1) notebook](simbolic_computations/MPRK_2_2_1_generalSystem.nb)
+* the computations of the $\Delta t$ bound and the Taylor expansion to understand which scheme has inconsistency on the **symmetric linear $2\times 2$ system** for some of the **MPRK(2,2,$\alpha$)** schemes is in the [MPRK(2,2,$\alpha$) notebook](simbolic_computations/MPRK_2_2_alpha.nb)
+* the Taylor expansion to understand if **MPRK(3,2)** has inconsistency on the **symmetric linear $2\times 2$ system** is in the [MPRK(3,2) notebook](simbolic_computations/MPRK_3_2.nb)
+* the Taylor expansion to understand if **mPDeC2** and **mPDeC3** have inconsistencies on the **symmetric linear $2\times 2$ system** is in the [MPDeC notebook](simbolic_computations/MPDeC.nb)
+
+### Supplementary material
+In folder [supplementary_material](supplementary_material/) there is the [file](supplementary_material/supplementary_material.pdf) with the supplementary results which are only described in the report 
+
+
+## Disclaimer
+
+Everything is provided as-is and without warranty. Use at your own risk!

simulations/Lagrange_polynomials.jl

+using SymPy
+
+
+function generate_poly(s, nodes,deg)
+#Generates Lagrangian polynomial in symbolic
+
+aux=[]
+phi=[]
+
+for k=1:deg+1
+    append!(aux, 1);
+    for j=1:deg+1
+        if j != k
+            aux[k]=aux[k].*(s-nodes[j]);
+        end
+    end
+end
+for j=1:deg+1
+    aux[j]= aux[j]/SymPy.N(subs(aux[j],s,nodes[j]));
+end
+return aux
+end
+
+function compute_theta(nodes)
+# Compute the integral of lagrange polynomials from 0 to interpolation point theta[r,m]= \int_{t^0}^{t^m} \phi_r
+s=Sym("s");
+deg=length(nodes)-1;
+poly=generate_poly(s,nodes,deg)
+
+theta=zeros(deg+1, deg+1)
+
+for r=1:deg+1
+    for m=1:deg+1
+        theta[r,m]=integrate(poly[r],(s,0,nodes[m]))
+    end
+end
+
+return theta
+
+end

simulations/mPDeC.jl

+using FastGaussQuadrature, Compat
+using LinearAlgebra
+
+"""Different models"""
+
+#include("Linear_model.jl")
+#include("Robertson.jl")
+#include("Nonlinear_model.jl")
+#include("Nonlinear_model2.jl")
+
+#include("Lagrange_polynomials.jl")
+#include("Error_module.jl")
+
+
+function lagrange_basis(X,t,i)
+    idxs = eachindex(X)
+    if t in X
+        if t == X[i]
+            return 1
+        else
+            return 0
+        end
+    end
+    return prod((t-X[j])/(X[i]-X[j]) for j in idxs if j != i)
+end
+
+
+function get_nodes(order,nodes_type)
+    if nodes_type=="equispaced"
+        nodes= range(0,stop=1, length=order)
+        w = 1/order*ones(order)
+    elseif nodes_type=="gaussLobatto"
+        nodes, w = gausslobatto(order)
+        w = map( x-> x*0.5, w)
+        nodes = nodes*0.5.+0.5
+    elseif nodes_type=="gaussLegendre"
+        nodes, w = gausslegendre(order)
+        w = map( x-> x*0.5, w)
+        nodes = nodes*0.5.+0.5
+    end
+    
+    return nodes, w
+end
+
+
+""" Computing the matrix of the values T_k^m = int_{t^0}^{t^m} varphi_k"""
+function compute_theta_DeC(order, nodes_type)
+    
+    nodes, w = get_nodes(order,nodes_type)
+    int_nodes, int_w = get_nodes(order,"gaussLobatto")
+    # generate theta coefficients 
+    theta = zeros(order,order)
+    for m in 1:order
+        nodes_m = int_nodes*(nodes[m])
+        w_m = int_w*(nodes[m])
+        for r in 1:order
+            theta[r,m] = sum([lagrange_basis(nodes,nodes_m[k],r)*w_m[k] for k in 1:order])
+        end
+    end
+    return theta
+    
+end
+
+
+#Dec_correction-
+#Func is the function of the ode
+#N_step -> Number of subtimesteps
+#Corr-> Number of Corrections
+#y_0 Initial condition
+# t_0 Start value
+# t_end End value
+# Maybe logical variable for checking for the different solvers or the distributaiton
+# of subtimesteps up to some order
+
+""" Dec_correction constant time step"""
+function dec_correction( prod_dest, rhs, tspan, y_0, M_sub::Int, K_corr::Int,
+    distribution)
+
+   N_time=length(tspan)
+   dim=length(y_0)
+   U=zeros(dim, N_time)
+   u_p=zeros(dim, M_sub+1)
+   u_a=zeros(dim, M_sub+1)
+   prod_p=zeros(dim,dim,M_sub+1)
+   dest_p=zeros(dim,dim,M_sub+1)
+   rhs_p =zeros(dim,M_sub+1)
+   Theta=compute_theta_DeC(M_sub+1,distribution)
+   U[:,1]=y_0
+   for it=2: N_time
+        delta_t=(tspan[it]-tspan[it-1])
+        for m=1:M_sub+1
+         u_a[:,m]=U[:,it-1]
+         u_p[:,m]=U[:,it-1]
+        end
+        for k=2:K_corr+1
+            u_p=copy(u_a)
+            for r=1:M_sub+1
+                prod_p[:,:,r], dest_p[:,:,r]=prod_dest(u_p[:,r])
+                rhs_p[:,r] = rhs(u_p[:,r])
+            end
+            for m=2:M_sub+1
+            u_a[:,m]=patanker_type_dec(prod_p, dest_p, rhs_p, delta_t, m, M_sub, Theta, u_p, dim)
+            end
+        end
+        U[:,it]=u_a[:,M_sub+1]
+    end
+    return tspan, U
+end
+
+
+
+"""prod_p and dest_p dim x dim
+delta_t is timestep length
+m_substep is the subtimesteps for which we are solving the system
+M_sub is the maximum number of subtimesteps (0,..., M)
+theta coefficient vector matrix  in  M_sub+1 x M_sub+1, first index is the index we are looping on,
+the second is the one we are solving for, theta_r,m= int_t0^tm phi_r
+u_p is the previous correction solution for all subtimesteps in dim x (M_sub+1) """
+function patanker_type_dec(prod_p, dest_p, rhs_p, delta_t,  m_substep::Int, M_sub, theta, u_p, dim)
+
+mass=Matrix{Float64}(I, dim, dim);
+#println(mass)
+RHS = u_p[:,1]
+for i=1:dim
+    for r=1: M_sub+1
+				RHS[i]=RHS[i] +delta_t*theta[r,m_substep]*rhs_p[i,r]
+        if theta[r,m_substep]>0
+            for j= 1: dim
+              mass[i,j]=mass[i,j]- delta_t*theta[r,m_substep]*(prod_p[i,j,r]/u_p[j,m_substep])
+              mass[i,i]=mass[i,i]+ delta_t*theta[r,m_substep]*(dest_p[i,j,r]/u_p[i,m_substep])
+            end
+        elseif theta[r,m_substep]<0
+            for j= 1: dim
+              mass[i,i]=mass[i,i]- delta_t*theta[r,m_substep]*(prod_p[i,j,r]/u_p[i,m_substep])
+              mass[i,j]=mass[i,j]+ delta_t*theta[r,m_substep]*(dest_p[i,j,r]/u_p[j,m_substep])
+            end
+        end
+    end
+end
+return mass\RHS
+end
+
+##dim=2;
+##M_sub=4;
+##u_p=rand(dim, M_sub+1 );
+##prod_p=dest_p= zeros(dim,dim,M_sub+1)
+##for r=1:M_sub+1
+##    prod_p[:,:,r],dest_p[:,:,r]=prod_dest(u_p[:,1])
+##end
+
+##delta_t=0.1
+##m_substep=2
+##theta=-1.0 .+2.0.*rand(M_sub+1,M_sub+1)
+
+##patanker_type_dec(prod_p,dest_p, delta_t, m_substep, M_sub, theta, u_p, dim)
+
+
+#####Convergence test
+###u0=[0.9;0.1]
+###dim = length(u0)
+###T_final = 1.75
+###n_ns = 5
+###Ns = 2 .^range(3,stop=2+n_ns)
+###dts = T_final./Ns
+###U_exact = exact_solution(u0,T_final)
+###plot()
+###for order = 2:6
+###  Us = zeros(dim, n_ns)
+###  for k=1:n_ns
+###     t,U=dec_correction(order,order+1,u0,0,T_final,Ns[k],true)
+###     Us[:,k]=U[:,end]
+###  end
+
+
+###  plot_convergence(dts,Us,U_exact,"dec $(order)", order)
+###end
+###plot!()
+
+
+
+"""
+Later Version for different quadrature rules, for later worth it to check!
+
+"""
+
+function quadrature_distribution_name_to_type(name::Symbol)
+    if name == :Periodic
+        #Going in the slope
+    elseif name == :Lobatto
+        gausslobatto
+    elseif name == :Chebyshev
+        gausschebyshev
+    else
+        error("Quadrature `$name` unknown.")
+    end
+end

simulations/mPRK.jl

+using LinearAlgebra
+
+function mPRK22(prod_dest,rhs, tspan, u0::AbstractArray{T}, α) where {T}
+    if α<1/2
+        error("Negative RK coefficients")
+    end
+    c=T[0,α,1] #c,1
+    A=T[[0 0]; [α 0]; [1-1/2/α 1/2/α]] #[A;b]
+    ex=T[[0 0 0];[1 0 0];[1-1/α 1/α 0]]
+    S = 3 # stages+1
+    dim=length(u0)
+    Nt=length(tspan)
+    U=zeros(T, dim,Nt)
+    y=zeros(T, S,dim)
+    p=zeros(T, S,dim,dim)
+    d=zeros(T, S,dim,dim)
+    rhs_p =zeros(T, S,dim)
+    dens=zeros(T, dim)
+    U[:,1]=u0
+    for it=2:Nt
+        dt=tspan[it]-tspan[it-1]
+        for k=1:S
+            rhs_p[k,:] = U[:,it-1]
+            for i=1:dim
+                dens[i]=1
+                for z=1:k-1
+                    dens[i]=dens[i]*y[z,i]^ex[k,z]
+                end
+            end
+            MM=Matrix{T}(I, dim, dim)
+            for z=1:k-1
+                rhs_p[k,:] = rhs_p[k,:] +dt*A[k,z]*rhs(y[z,:])
+                for i=1:dim
+                    for j=1:dim
+                        MM[i,i]=MM[i,i]+A[k,z]*dt*d[z,i,j]/dens[i]
+                        MM[i,j]=MM[i,j]-A[k,z]*dt*p[z,i,j]/dens[j]
+                    end
+                end
+            end
+            y[k,:]=MM\rhs_p[k,:]
+            p[k,:,:],d[k,:,:]=prod_dest(y[k,:])
+        end
+        U[:,it]=y[S,:]
+    end
+    
+    return tspan, U
+end
+
+
+function mPRK32(prod_dest,rhs, tspan, u0)
+    c=[0,1,0.5,1] #c,1
+    A=[[0 0 0]; [1 0 0]; [0.25 0.25 0]; [1/6 1/6 2/3]] #[A;b]
+    S = 4 # stages+1
+    dim=length(u0)
+    Nt=length(tspan)
+    U=zeros(dim,Nt)
+    y=zeros(S,dim)
+    p=zeros(S,dim,dim)
+    d=zeros(S,dim,dim)
+    rhs_p =zeros(S,dim)
+    dens=zeros(dim)
+    U[:,1]=u0
+    for it=2:Nt
+        dt=tspan[it]-tspan[it-1]
+        for k=1:S
+            rhs_p[k,:] = U[:,it-1]
+            for i=1:dim
+                if k==1
+                    dens[i]=1
+                else
+                    dens[i]=y[minimum([k-1,2]),i]
+                end
+            end
+            MM=Matrix{Float64}(I, dim, dim)
+            for z=1:k-1
+                rhs_p[k,:] = rhs_p[k,:] +dt*A[k,z]*rhs(y[z,:])
+                for i=1:dim
+                    for j=1:dim
+                        MM[i,i]=MM[i,i]+A[k,z]*dt*d[z,i,j]/dens[i]
+                        MM[i,j]=MM[i,j]-A[k,z]*dt*p[z,i,j]/dens[j]
+                    end
+                end
+            end
+            y[k,:]=MM\rhs_p[k,:]
+            p[k,:,:],d[k,:,:]=prod_dest(y[k,:])
+        end
+        U[:,it]=y[S,:]
+    end
+    
+    return tspan, U
+end
+
+
+function mPRKSO2(prod_dest,rhs,tspan, u0,α,β)
+    if α<0 || α>1 || β<0 || α*β+1/2/β>1 
+        error("Negative RK coefficients")
+    end
+    c=[0,α,1] #c,1
+    A=[[0 0]; [1 0]; [1-α α]] #[A;b]
+    b = [[0 0]; [β 0]; [1-1/2/β-β*α 1/2/β]]
+    s=(1-α*β+α*β^2)/β/(1-α*β)
+    ex=[[0 0 0];[1 0 0];[1-s s 0]]
+    S = 3 # stages+1
+    dim=length(u0)
+    Nt=length(tspan)
+    U=zeros(dim,Nt)
+    y=zeros(S,dim)
+    p=zeros(S,dim,dim)
+    d=zeros(S,dim,dim)
+    rhs_p =zeros(S,dim)
+    dens=zeros(dim)
+    U[:,1]=u0
+    for it=2:Nt
+        dt=tspan[it]-tspan[it-1]
+        y[1,:]=U[:,it-1]
+				p[1,:,:],d[1,:,:]=prod_dest(y[1,:])
+        rhs_p[1,:] = zeros(dim)
+        for k=2:S
+        rhs_p[k,:] = zeros(dim)
+            for i=1:dim
+                dens[i]=1
+                for z=1:k-1
+                    dens[i]=dens[i]*y[z,i]^ex[k,z]
+                end
+            end
+            MM=Matrix{Float64}(I, dim, dim)
+            for z=1:k-1
+                rhs_p[k,:] = rhs_p[k,:] +dt*b[k,z]*rhs(y[z,:])
+				        for i=1:dim
+				            for j=1:dim
+                        MM[i,i]=MM[i,i]+b[k,z]*dt*d[z,i,j]/dens[i]
+                        MM[i,j]=MM[i,j]-b[k,z]*dt*p[z,i,j]/dens[j]
+                    end
+                end
+            end
+            RHS=A[:,:]*y[1:S-1,:]
+            RHS[k,:] = RHS[k,:]+rhs_p[k,:]
+            y[k,:]=MM\RHS[k,:]
+            p[k,:,:],d[k,:,:]=prod_dest(y[k,:])
+        end
+        U[:,it]=y[S,:]
+    end
+    
+    return tspan, U
+end
+
+
+function mPRK3(prod_dest,rhs,tspan, u0,α,β)
+    α0=1/6*(3+(3-2*sqrt(2))^(1/3)+(3+2*sqrt(2))^(1/3))
+    if α<1/3  || β<0
+        error("Negative RK coefficients")
+    elseif α<2/3 && (β<2/3 || β >3*α*(1-α))
+        error("Negative RK coefficients")
+    elseif α>=2/3 && α<α0 && (β>2/3 || β <3*α*(1-α))
+        error("Negative RK coefficients")
+    elseif α>=α0 && (β>2/3 || β <(3*α-2)/(6*α-3))
+        error("Negative RK coefficients")
+    end
+    
+    c=[0,α,β, 1]
+    A=[ [0 0 0]; [α 0 0]; [(3*α*β*(1-α)-β^2)/α/(2-3*α)  β*(β-α)/α/(2-3*α) 0]; [1+(2-3*(α+β))/6/α/β  (3*β-2)/(6*α*(β-α)) (2-3*α)/(6*β*(β-α))]]
+    if maximum(A)<0 || maximum(c)<0
+        error("Negative RK coefficients")
+    end
+
+    S = 4 # stages+1
+
+    pp=3*A[2,1]*(A[3,1]+A[3,2])*A[S,3]
+    q=A[2,1]
+    beta=zeros(2)
+    beta[2]=1/2/A[2,1]
+    beta[1]=1-beta[2]
+
+    dim=length(u0)
+    Nt=length(tspan)
+    U=zeros(dim,Nt)
+    y=zeros(S,dim)
+    p=zeros(S,dim,dim)
+    d=zeros(S,dim,dim)
+    rhs_p =zeros(S,dim)
+    dens=zeros(dim)
+    U[:,1]=u0
+    for it=2:Nt
+        dt=tspan[it]-tspan[it-1]
+        # u1
+        y[1,:]=U[:,it-1]
+        p[1,:,:],d[1,:,:]=prod_dest(y[1,:])
+        rhs_p[1,:] = rhs(y[1,:])
+
+        # u2
+        k=2
+        rhs_p[k,:] = zeros(dim)
+        MM=Matrix{Float64}(I, dim, dim)
+        for z=1:k-1
+            rhs_p[k,:] = rhs_p[k,:] +dt*A[k,z]*rhs(y[z,:])
+            for i=1:dim
+                for j=1:dim
+                    MM[i,i]=MM[i,i]+A[k,z]*dt*d[z,i,j]/y[1,i]
+                    MM[i,j]=MM[i,j]-A[k,z]*dt*p[z,i,j]/y[1,j]
+                end
+            end
+        end
+        RHS= y[1,:]+rhs_p[k,:]
+        y[k,:]=MM\RHS
+        p[k,:,:],d[k,:,:]=prod_dest(y[k,:])
+        rhs_p[k,:] = rhs(y[k,:])
+
+        # u3
+        k=3
+        rhs_p[k,:] = zeros(dim)
+        rho=zeros(dim)
+
+        for i=1:dim
+            rho[i]=y[2,i]^(1/pp)*y[1,i]^(1-1/pp)
+        end
+
+        MM=Matrix{Float64}(I, dim, dim)
+        for z=1:k-1
+            rhs_p[k,:] = rhs_p[k,:] +dt*A[k,z]*rhs(y[z,:])
+            for i=1:dim
+                for j=1:dim
+                    MM[i,i]=MM[i,i]+A[k,z]*dt*d[z,i,j]/rho[i]
+                    MM[i,j]=MM[i,j]-A[k,z]*dt*p[z,i,j]/rho[j]
+                end
+            end
+        end
+        RHS = y[1,:]+rhs_p[k,:]
+        y[k,:]=MM\RHS
+        p[k,:,:],d[k,:,:]=prod_dest(y[k,:])
+        rhs_p[k,:] = rhs(y[k,:])
+
+        # u4
+        k=4
+                        
+        mu=zeros(dim)
+        as = zeros(dim)
+
+        for i=1:dim
+            mu[i]=U[i,it-1].^(1-1/q).*y[2,i].^(1/q)
+        end
+        
+        rrr=zeros(dim)
+        MM=Matrix{Float64}(I, dim, dim)
+        rrr = y[1,:] + dt*beta[1]*rhs(y[1,:])+dt*beta[2]*rhs(y[2,:])
+        for i = 1:dim
+            for j=1:dim
+                 MM[i,j] = MM[i,j] -dt*(beta[1]*p[1,i,j]+beta[2]*p[2,i,j])/mu[j]
+                 MM[i,i] = MM[i,i] +dt*(beta[1]*d[1,i,j]+beta[2]*d[2,i,j])/mu[i]
+            end
+        end
+        σ=MM\rrr
+
+        rhs_p[k,:] = zeros(dim)
+        MM=Matrix{Float64}(I, dim, dim)
+        for z=1:k-1
+            rhs_p[k,:] = rhs_p[k,:] +dt*A[k,z]*rhs(y[z,:])
+            for i=1:dim
+                for j=1:dim
+                    MM[i,i]=MM[i,i]+A[k,z]*dt*d[z,i,j]/σ[i]
+                    MM[i,j]=MM[i,j]-A[k,z]*dt*p[z,i,j]/σ[j]
+                end
+            end
+        end
+        RHS = y[1,:]+rhs_p[k,:]
+        y[k,:]=MM\RHS
+        p[k,:,:],d[k,:,:]=prod_dest(y[k,:])
+        U[:,it]=y[S,:]
+    end
+    return tspan, U
+end
+
+
+function mPRKSO3(prod_dest,rhs,tspan, u0)
+    A=[[0 0 0]; [1 0 0]; [9.2600312554031827E-01 7.3996874459681783E-02 0]; 
+       [7.0439040373427619E-01 2.0662904223744017E-10 2.9560959605909481E-01]] #[A;b]
+    b = [[0 0 0 0]; [4.7620819268131703E-01 0 0 0]; [7.7545442722396801E-02 5.9197500149679749E-01 0 0];
+         [2.0044747790361456E-01 6.8214380786704851E-10 5.9121918658514827E-01 0]]
+
+    n1=2.569046025732011E-01; n2 = 7.430953974267989E-01; z = 6.288938077828750E-01;	
+    η1 = 3.777285888379173E-02; η2 = 1/3; η3 = 1.868649805549811E-01; η4=2.224876040351123; s = 5.721964308755304
+
+    S = 4 # stages+1
+    dim=length(u0)
+    Nt=length(tspan)
+    U=zeros(dim,Nt)
+    y=zeros(S,dim)
+    p=zeros(S,dim,dim)