ODE15s with non-constant Jacobian

Hi everyone,

i want to solve an ODE by using the Jacobian, that is not constant:

options = odeset( 'Jacobian' , @(t,x)J(u_int, deltaZ, Nz, cFeed), 'RelTol',1e-8, 'AbsTol',1e-8);

[t,V] = ode15s(@(t,v)odeTest(t,v), tspan, v0, options);

function dfdx = J(u_int, deltaZ, Nz, cFeed)

x = sym('x' , [1 Nz+4]);

f1 = (-u_int .* (x(1)^2 - cFeed) ./ deltaZ );

f2 = (-u_int .* (x(2:Nz+4) - x(1:Nz+3)) ./ deltaZ );

f = [f1, f2];

J = jacobian(f,x);

dfdx = (J);

end

%% ODE

function dfdt = odeTest(~, vi)

dfdt = zeros(length(vi),1);

%% Variables

u_int = 0.10;

H = 0.20;

Nz = 100;

deltaZ = H/(Nz+4);

cFeed = 10;

x = vi(1:Nz+4);

%% DGL

dfdt(1) = (-u_int .* (x(1)^2 - cFeed) ./ deltaZ );

dfdt(2: Nz+4) = (-u_int .* (x(2:Nz+4) - x(1:Nz+3)) ./ deltaZ );

end

The problem I am facing is following error:

Conversion to logical from sym is not possible.

Error in ode15s (line 405)

if absh * rh > 1

Error in test_jacobi_upwind (line 39)

[t,V] = ode15s(@(t,v)odeTest(t,v), tspan, v0, options);

I can make it work, if the Jacobian is a constant and does not depend on x(i) - but as soon as J is depending on x(i), I am getting errors.

Any help is appreciated.

Best,

M

Answer 1

[Torsten]

JacFun = J(u_int, deltaZ, Nz, cFeed);
options = odeset( 'Jacobian' , @(t,x)JacFun(t, x, u_int, deltaZ, Nz, cFeed), 'RelTol',1e-8, 'AbsTol',1e-8);
[t,V] = ode15s(@(t,x)odeTest(t,x,u_int, deltaZ, Nz, cFeed), tspan, v0, options);

function dfdx = J(u_int, deltaZ, Nz, cFeed)
  x = sym('x' , [1 Nz+4]);
  f1 = (-u_int .*  (x(1)^2      - cFeed) ./ deltaZ );
  f2 = (-u_int .*  (x(2:Nz+4)   - x(1:Nz+3)) ./ deltaZ );
  f = [f1, f2];
  J = jacobian(f,x);
  dfdx = matlabFunction(J,'Vars',[t, x, u_int, deltaZ, Nz, cFeed]);
end

function dfdt = odeTest(t, vi, u_int, deltaZ, Nz, cFeed)       
  dfdt = zeros(length(vi),1);   
  %% Variables
  u_int  = 0.10;
  H      = 0.20;
  Nz     = 100;
  deltaZ = H/(Nz+4);
  cFeed  = 10;

  x   = vi(1:Nz+4);
  %% DGL
  dfdt(1)              = (-u_int .*  (x(1)^2   - cFeed) ./ deltaZ );
  dfdt(2: Nz+4)        = (-u_int .*  (x(2:Nz+4)   - x(1:Nz+3)) ./ deltaZ );
end

[Berry]

clear

clc

%% Settings

Nz = 100;

time = 100;

total_num = Nz + 4;

v0 = zeros(total_num,1);

tspan = (0:0.01:time);

%%Parameter

u_int = 0.10;

H = 0.20;

deltaZ = H/(Nz+4);

cFeed = 10;

%% Torsten

JacFun = J(u_int, deltaZ, Nz, cFeed);

options = odeset( 'Jacobian' , @(t,x)JacFun(t, x, u_int, deltaZ, Nz, cFeed), 'RelTol',1e-8, 'AbsTol',1e-8);

[t,V] = ode15s(@(t,x)odeTest(t,x,u_int, deltaZ, Nz, cFeed), tspan, v0, options);

function dfdx = J(u_int, deltaZ, Nz, cFeed)

x = sym('x' , [1 Nz+4]);

f1 = (-u_int .* (x(1)^2 - cFeed) ./ deltaZ );

f2 = (-u_int .* (x(2:Nz+4) - x(1:Nz+3)) ./ deltaZ );

f = [f1, f2];

J = jacobian(f,x);

dfdx = matlabFunction(J,'Vars',[t, x, u_int, deltaZ, Nz, cFeed]);

end

function dfdt = odeTest(t, vi, u_int, deltaZ, Nz, cFeed)

dfdt = zeros(length(vi),1);

%% Variables

u_int = 0.10;

H = 0.20;

Nz = 100;

deltaZ = H/(Nz+4);

cFeed = 10;

x = vi(1:Nz+4);

%% DGL

dfdt(1) = (-u_int .* (x(1)^2 - cFeed) ./ deltaZ );

dfdt(2: Nz+4) = (-u_int .* (x(2:Nz+4) - x(1:Nz+3)) ./ deltaZ );

end

Dear @Torsten,

thank you! Unfortunately, it is still not working:

Unrecognized function or variable 't'.

Error in test_jacobi_upwind_MATHWORKS>J (line 31)

dfdx = matlabFunction(J,'Vars',[t, x, u_int, deltaZ, Nz, cFeed]);

Error in test_jacobi_upwind_MATHWORKS (line 21)

JacFun = J(u_int, deltaZ, Nz, cFeed);

I tried to delete "t" everywhere, but that is not the solution.

[Torsten]

clear
clc
%% Settings
Nz     = 100;
time = 100;
total_num = Nz + 4;
v0 = zeros(total_num,1);
tspan = (0:0.01:time);
%%Parameter
u_int  = 0.10;
H      = 0.20;
deltaZ = H/(Nz+4);
cFeed  = 10;
%% Torsten
JacFun = J(u_int, deltaZ, Nz, cFeed);
options = odeset( 'Jacobian' , JacFun, 'RelTol',1e-8, 'AbsTol',1e-8);
[t,V] = ode15s(@(t,x)odeTest(t,x,u_int, deltaZ, Nz, cFeed), tspan, v0, options);
function dfdx = J(u_int, deltaZ, Nz, cFeed)
  syms t
  x = sym('x' , [1 Nz+4]);
  f1 = (-u_int .*  (x(1)^2      - cFeed) ./ deltaZ );
  f2 = (-u_int .*  (x(2:Nz+4)   - x(1:Nz+3)) ./ deltaZ );
  f = [f1, f2];
  J = jacobian(f,x);
  dfdx = matlabFunction(J,'Vars',{t, [x.']});
end
function dfdt = odeTest(t, vi, u_int, deltaZ, Nz, cFeed)       
  dfdt = zeros(length(vi),1);   
  %% Variables
  u_int  = 0.10;
  H      = 0.20;
  Nz     = 100;
  deltaZ = H/(Nz+4);
  cFeed  = 10;
  x   = vi(1:Nz+4);
  %% DGL
  dfdt(1)              = (-u_int .*  (x(1)^2   - cFeed) ./ deltaZ );
  dfdt(2: Nz+4)        = (-u_int .*  (x(2:Nz+4)   - x(1:Nz+3)) ./ deltaZ );
end

[Berry]

Thank you @Torsten! Now it worked.

So now I wanted to use this for the WENO scheme instead of the upwind scheme (see this post: ODE solver with WENO scheme (weighted essential non-oscillatory) - (mathworks.com))

clear

clc

%% Settings

Nz = 20;

time = 100;

total_num = Nz + 4;

v0 = zeros(total_num,1);

tspan = (0:0.01:time);

%% Variables

cFeed = 10;

u_int = 0.10;

H = 0.20;

deltaZ = H/(Nz+4);

epsilon = 1e-12;

%% Torsten

JacFunc = J(u_int, deltaZ, Nz, epsilon, cFeed);

options = odeset( 'Jacobian' , @(t,x)JacFunc(u_int, deltaZ, Nz, epsilon,cFeed), 'RelTol',1e-8, 'AbsTol',1e-8);

%options = odeset( 'RelTol',1e-8, 'AbsTol',1e-8);

[t,V] = ode15s(@(t,x)odeTest(t,x, u_int, deltaZ, Nz,epsilon, cFeed), tspan, v0, options);

dbk = [t,V];

%% Plotting

figure('Name', 'Upwind_DBK')

plot(dbk(:,1), dbk(:,(Nz+4)));

%% Jacobian

function dfdx = J(u_int, deltaZ, Nz, epsilon, cFeed)

syms t

x = sym('x' , [1 Nz+4]);

f1 = -u_int * (x(1) - cFeed) ./ deltaZ;

f2 = -u_int * (...

( (...

(2/3 /( epsilon+(x(3) - x(2) ).^2 ).^2 ) ...

./ ( (2/3 /( epsilon+(x(3) - x(2) ).^2 ).^2 ) ...

+ (1/3 /( epsilon+(x(2) - x(1) ).^2 ).^2 ) ...

)...

) ...

* ( 0.5 * x(2) + 0.5 * x(3)) ...

+ (...

(1/3 /( epsilon+(x(2) - x(1) ).^2 ).^2 ) ...

./ ( (2/3 /( epsilon+(x(3) - x(2) ).^2 ).^2 ) ...

+ (1/3 /( epsilon+(x(2) - x(1) ).^2 ).^2 ) ...

)...

) ...

* (-1/2 * x(1) + 3/2 * x(2)) ...

) ...

- ...

( (...

(2/3 /( epsilon+(x(2) - x(1) ).^2 ).^2 ) ...

./ ( (2/3 /( epsilon+(x(2) - x(1) ).^2 ).^2 ) ...

+ (1/3 /( epsilon+(x(1) - cFeed ).^2 ).^2 ) ...

)...)

) ...

* ( 0.5 * x(1) + 0.5 * x(2)) ...

+ (...

(1/3 /( epsilon+(x(1) - cFeed ).^2 ).^2 ) ...

./ ( (2/3 /( epsilon+(x(2) - x(1) ).^2 ).^2 ) ...

+ (1/3 /( epsilon+(x(1) - cFeed ).^2 ).^2 ) ...

)...

) ...

* (-1/2 * cFeed + 3/2 * x(1)) ...

) ...

)...

./ deltaZ;

f3 = -u_int *( ...

( ( ...

(0.3 ./( epsilon + ((13/12.*( x(3)- 2.*x(4) + x(5)).^2)+ ( 1/4.*( 3.*x(2)- 4.*x(4) + x(5) ) .^2))).^2) ...

./ ( (0.3 ./( epsilon + ((13/12.*( x(3)- 2.*x(4) + x(5)).^2)+ ( 1/4.*( 3.*x(2)- 4.*x(4) + x(5) ) .^2))).^2) ...

+ (0.6 ./( epsilon + ((13/12.*( x(2)- 2.*x(3) + x(4)).^2)+ ( 1/4.*( x(2)- x(4) ) .^2))).^2) ...

+ (0.1 ./( epsilon + ((13/12.*( x(1)- 2.*x(2) + x(3)).^2)+ ( 1/4.*( x(1)- 4.*x(2) + 3*x(3) ) .^2))).^2) ...

) ...

) ...

.* ( 1/3 .* x(3) + 5/6 .*x(4) - 1/6 .*x(5)) ...

+ (...

(0.6 ./( epsilon + ((13/12.*( x(2)- 2.*x(3) + x(4)).^2)+ ( 1/4.*( x(2)- x(4) ) .^2))).^2) ...

./ ( (0.3 ./( epsilon + ((13/12.*( x(3)- 2.*x(4) + x(5)).^2)+ ( 1/4.*( 3.*x(2)- 4.*x(4) + x(5) ) .^2))).^2) ...

+ (0.6 ./( epsilon + ((13/12.*( x(2)- 2.*x(3) + x(4)).^2)+ ( 1/4.*( x(2)- x(4) ) .^2))).^2) ...

+ (0.1 ./( epsilon + ((13/12.*( x(1)- 2.*x(2) + x(3)).^2)+ ( 1/4.*( x(1)- 4.*x(2) + 3*x(3) ) .^2))).^2) ...

)...

) ...

.* (- 1/6 .* x(2) + 5/6 .*x(3) + 1/3 .*x(4)) ...

+ (...

(0.1 ./( epsilon + ((13/12.*( x(1)- 2.*x(2) + x(3)).^2)+ ( 1/4.*( x(1)- 4.*x(2) + 3*x(3) ) .^2))).^2) ...

./ ( (0.3 ./( epsilon + ((13/12.*( x(3)- 2.*x(4) + x(5)).^2)+ ( 1/4.*( 3.*x(2)- 4.*x(4) + x(5) ) .^2))).^2) ...

+ (0.6 ./( epsilon + ((13/12.*( x(2)- 2.*x(3) + x(4)).^2)+ ( 1/4.*( x(2)- x(4) ) .^2))).^2) ...

+ (0.1 ./( epsilon + ((13/12.*( x(1)- 2.*x(2) + x(3)).^2)+ ( 1/4.*( x(1)- 4.*x(2) + 3*x(3) ) .^2))).^2) ...

)...

) ...

.* ( 1/3 .* x(1) - 7/6 .*x(2) + 11/6 .*x(3)) ...

)...

-...

( (...

(2/3 /( epsilon+(x(3) - x(2) ).^2 ).^2 ) ...

./ ( (2/3 /( epsilon+(x(3) - x(2) ).^2 ).^2 ) ...

+ (1/3 /( epsilon+(x(2) - x(1) ).^2 ).^2 ) ...

)...

) ...

* ( 0.5 * x(2) + 0.5 * x(3)) ...

+ (...

(1/3 /( epsilon+(x(2) - x(1) ).^2 ).^2 ) ...

./ ( (2/3 /( epsilon+(x(3) - x(2) ).^2 ).^2 ) ...

+ (1/3 /( epsilon+(x(2) - x(1) ).^2 ).^2 ) ...

)...

) ...

* (-1/2 * x(1) + 3/2 * x(2)) ...

)...

) ./ deltaZ;

f4 = -u_int *( ...

( ( ...

(0.3 ./( epsilon + ((13/12.*( x(4:Nz+2)- 2.*x(5:Nz+3) + x(6:Nz+4)).^2)+ ( 1/4.*( 3.*x(3:Nz+1)- 4.*x(5:Nz+3) + x(6:Nz+4) ).^2))).^2) ...

./ ( (0.3 ./( epsilon + ((13/12.*( x(4:Nz+2)- 2.*x(5:Nz+3) + x(6:Nz+4)).^2)+ ( 1/4.*( 3.*x(3:Nz+1)- 4.*x(5:Nz+3) + x(6:Nz+4) ).^2))).^2) ...

+ (0.6 ./( epsilon + ((13/12.*( x(3:Nz+1)- 2.*x(4:Nz+2) + x(5:Nz+3)).^2)+ ( 1/4.*( x(3:Nz+1)- x(5:Nz+3) ).^2))).^2) ...

+ (0.1 ./( epsilon + ((13/12.*( x(2:Nz) - 2.*x(3:Nz+1) + x(4:Nz+2)).^2)+ ( 1/4.*( x(2:Nz) - 4.*x(3:Nz+1) + 3*x(4:Nz+2) ).^2))).^2) ...

) ...

) ...

.* ( 1/3 .* x(4:Nz+2) + 5/6 .*x(5:Nz+3) - 1/6 .*x(6:Nz+4)) ...

+ (...

(0.6 ./( epsilon + ((13/12.*( x(3:Nz+1)- 2.*x(4:Nz+2) + x(5:Nz+3)).^2)+ ( 1/4.*( x(3:Nz+1)- x(5:Nz+3) ).^2))).^2) ...

./ ( (0.3 ./( epsilon + ((13/12.*( x(4:Nz+2)- 2.*x(5:Nz+3) + x(6:Nz+4)).^2)+ ( 1/4.*( 3.*x(3:Nz+1)- 4.*x(5:Nz+3) + x(6:Nz+4) ).^2))).^2) ...

+ (0.6 ./( epsilon + ((13/12.*( x(3:Nz+1)- 2.*x(4:Nz+2) + x(5:Nz+3)).^2)+ ( 1/4.*( x(3:Nz+1)- x(5:Nz+3) ).^2))).^2) ...

+ (0.1 ./( epsilon + ((13/12.*( x(2:Nz) - 2.*x(3:Nz+1) + x(4:Nz+2)).^2)+ ( 1/4.*( x(2:Nz) - 4.*x(3:Nz+1) + 3*x(4:Nz+2) ).^2))).^2) ...

)...

) ...

.* (- 1/6 .* x(3:Nz+1) + 5/6 .*x(4:Nz+2) + 1/3 .*x(5:Nz+3)) ...

+ (...

(0.1 ./( epsilon + ((13/12.*( x(2:Nz) - 2.*x(3:Nz+1) + x(4:Nz+2)).^2)+ ( 1/4.*( x(2:Nz) - 4.*x(3:Nz+1) + 3*x(4:Nz+2) ).^2))).^2) ...

./ ( (0.3 ./( epsilon + ((13/12.*( x(4:Nz+2)- 2.*x(5:Nz+3) + x(6:Nz+4)).^2)+ ( 1/4.*( 3.*x(3:Nz+1)- 4.*x(5:Nz+3) + x(6:Nz+4) ).^2))).^2) ...

+ (0.6 ./( epsilon + ((13/12.*( x(3:Nz+1)- 2.*x(4:Nz+2) + x(5:Nz+3)).^2)+ ( 1/4.*( x(3:Nz+1)- x(5:Nz+3) ).^2))).^2) ...

+ (0.1 ./( epsilon + ((13/12.*( x(2:Nz) - 2.*x(3:Nz+1) + x(4:Nz+2)).^2)+ ( 1/4.*( x(2:Nz) - 4.*x(3:Nz+1) + 3*x(4:Nz+2) ).^2))).^2) ...

)...

) ...

.* ( 1/3 .* x(2:Nz) - 7/6 .*x(3:Nz+1) + 11/6 .*x(4:Nz+2)) ...

)...

-...

( ( ...

(0.3 ./( epsilon + ((13/12.*( x(3:Nz+1)- 2.*x(4:Nz+2) + x(5:Nz+3)).^2)+ ( 1/4.*( 3.*x(2:Nz) - 4.*x(4:Nz+2) + x(5:Nz+3) ).^2))).^2) ...

./ ( (0.3 ./( epsilon + ((13/12.*( x(3:Nz+1)- 2.*x(4:Nz+2) + x(5:Nz+3)).^2)+ ( 1/4.*( 3.*x(2:Nz) - 4.*x(4:Nz+2) + x(5:Nz+3) ).^2))).^2) ...

+ (0.6 ./( epsilon + ((13/12.*( x(2:Nz) - 2.*x(3:Nz+1) + x(4:Nz+2)).^2)+ ( 1/4.*( x(2:Nz) - x(4:Nz+2) ).^2))).^2) ...

+ (0.1 ./( epsilon + ((13/12.*( x(1:Nz-1)- 2.*x(2:Nz) + x(3:Nz+1)).^2)+ ( 1/4.*( x(1:Nz-1)- 4.*x(2:Nz) + 3*x(3:Nz+1) ).^2))).^2) ...

) ...

) ...

.* ( 1/3 .* x(3:Nz+1) + 5/6 .*x(4:Nz+2) - 1/6 .*x(5:Nz+3)) ...

+ (...

(0.6 ./( epsilon + ((13/12.*( x(2:Nz) - 2.*x(3:Nz+1) + x(4:Nz+2)).^2)+ ( 1/4.*( x(2:Nz) - x(4:Nz+2) ).^2))).^2) ...

./ ( (0.3 ./( epsilon + ((13/12.*( x(3:Nz+1)- 2.*x(4:Nz+2) + x(5:Nz+3)).^2)+ ( 1/4.*( 3.*x(2:Nz) - 4.*x(4:Nz+2) + x(5:Nz+3) ).^2))).^2) ...

+ (0.6 ./( epsilon + ((13/12.*( x(2:Nz) - 2.*x(3:Nz+1) + x(4:Nz+2)).^2)+ ( 1/4.*( x(2:Nz) - x(4:Nz+2) ).^2))).^2) ...

+ (0.1 ./( epsilon + ((13/12.*( x(1:Nz-1)- 2.*x(2:Nz) + x(3:Nz+1)).^2)+ ( 1/4.*( x(1:Nz-1)- 4.*x(2:Nz) + 3*x(3:Nz+1) ).^2))).^2) ...

)...

) ...

.* (- 1/6 .* x(2:Nz) + 5/6 .*x(3:Nz+1) + 1/3 .*x(4:Nz+2)) ...

+ (...

(0.1 ./( epsilon + ((13/12.*( x(1:Nz-1) - 2.*x(2:Nz) + x(3:Nz+1)).^2)+ ( 1/4.*( x(1:Nz-1)- 4.*x(2:Nz) + 3*x(3:Nz+1) ).^2))).^2) ...

./ ( (0.3 ./( epsilon + ((13/12.*( x(3:Nz+1)- 2.*x(4:Nz+2) + x(5:Nz+3)).^2)+ ( 1/4.*( 3.*x(2:Nz) - 4.*x(4:Nz+2) + x(5:Nz+3) ).^2))).^2) ...

+ (0.6 ./( epsilon + ((13/12.*( x(2:Nz) - 2.*x(3:Nz+1) + x(4:Nz+2)).^2)+ ( 1/4.*( x(2:Nz) - x(4:Nz+2) ).^2))).^2) ...

+ (0.1 ./( epsilon + ((13/12.*( x(1:Nz-1)- 2.*x(2:Nz) + x(3:Nz+1)).^2)+ ( 1/4.*( x(1:Nz-1)- 4.*x(2:Nz) + 3*x(3:Nz+1) ).^2))).^2) ...

)...

) ...

.* ( 1/3 .* x(1:Nz-1) - 7/6 .*x(2:Nz) + 11/6 .*x(3:Nz+1)) ...

)...

)...

./ deltaZ;

f5 = -u_int * (...

( (...

(2/3 /( epsilon+(x(Nz+4) - x(Nz+3) ).^2 ).^2 ) ...

./ ( (2/3 /( epsilon+(x(Nz+4) - x(Nz+3) ).^2 ).^2 ) ...

+ (1/3 /( epsilon+(x(Nz+3) - x(Nz+2) ).^2 ).^2 ) ...

)...

) ...

* ( 0.5 * x(Nz+3) + 0.5 * x(Nz+4)) ...

+ (...

(1/3 /( epsilon+(x(Nz+3) - x(Nz+2) ).^2 ).^2 ) ...

./ ( (2/3 /( epsilon+(x(Nz+4) - x(Nz+3) ).^2 ).^2 ) ...

+ (1/3 /( epsilon+(x(Nz+3) - x(Nz+2) ).^2 ).^2 ) ...

)...

) ...

* (-1/2 * x(Nz+2) + 3/2 * x(Nz+3)) ...

) ...

- ...

( (...

(2/3 /( epsilon+(x(Nz+3) - x(Nz+2) ).^2 ).^2 ) ...

./ ( (2/3 /( epsilon+(x(Nz+3) - x(Nz+2) ).^2 ).^2 ) ...

+ (1/3 /( epsilon+(x(Nz+2) - x(Nz+1) ).^2 ).^2 ) ...

)...

) ...

* ( 0.5 * x(Nz+2) + 0.5 * x(Nz+3)) ...

+ (...

(1/3 /( epsilon+(x(Nz+2) - x(Nz+1) ).^2 ).^2 ) ...

./ ( (2/3 /( epsilon+(x(Nz+3) - x(Nz+2) ).^2 ).^2 ) ...

+ (1/3 /( epsilon+(x(Nz+2) - x(Nz+1) ).^2 ).^2 ) ...

)...

) ...

* (-1/2 * x(Nz+1) + 3/2 * x(Nz+2)) ...

)...

)./ deltaZ;

f6 = -u_int * (x(Nz+4) - x(Nz+3)) ./ deltaZ;

f = [f1, f2, f3, f4, f5, f6];

J = jacobian(f,x);

dfdx = matlabFunction(J,'Vars',{t, [x.']});

end

%% DGLs

function dfdt = odeTest(t, vi, u_int, deltaZ, Nz, epsilon, cFeed)

dfdt = zeros(length(vi),1);

%% Variables

Nz = 20;

epsilon = 1e-12;

u_int = 0.10;

H = 0.20;

deltaZ = H/(Nz+4);

cFeed = 10;

x = vi(1:Nz+4);

%% DGL

%% Upwind - start

dfdt(1) = -u_int * (x(1) - cFeed) ./ deltaZ;

%% WENO k=2 start

dfdt(2) = -u_int * (...

( (...

(2/3 /( epsilon+(x(3) - x(2) ).^2 ).^2 ) ...

./ ( (2/3 /( epsilon+(x(3) - x(2) ).^2 ).^2 ) ...

+ (1/3 /( epsilon+(x(2) - x(1) ).^2 ).^2 ) ...

)...

) ...

* ( 0.5 * x(2) + 0.5 * x(3)) ...

+ (...

(1/3 /( epsilon+(x(2) - x(1) ).^2 ).^2 ) ...

./ ( (2/3 /( epsilon+(x(3) - x(2) ).^2 ).^2 ) ...

+ (1/3 /( epsilon+(x(2) - x(1) ).^2 ).^2 ) ...

)...

) ...

* (-1/2 * x(1) + 3/2 * x(2)) ...

) ...

- ...

( (...

(2/3 /( epsilon+(x(2) - x(1) ).^2 ).^2 ) ...

./ ( (2/3 /( epsilon+(x(2) - x(1) ).^2 ).^2 ) ...

+ (1/3 /( epsilon+(x(1) - cFeed ).^2 ).^2 ) ...

)...

) ...

* ( 0.5 * x(1) + 0.5 * x(2)) ...

+ (...

(1/3 /( epsilon+(x(1) - cFeed ).^2 ).^2 ) ...

./ ( (2/3 /( epsilon+(x(2) - x(1) ).^2 ).^2 ) ...

+ (1/3 /( epsilon+(x(1) - cFeed ).^2 ).^2 ) ...

)...

) ...

* (-1/2 * cFeed + 3/2 * x(1)) ...

) ...

)./ deltaZ;

%% 3

dfdt(3) = -u_int *( ...

( ( ...

(0.3 ./( epsilon + ((13/12.*( x(3)- 2.*x(4) + x(5)).^2)+ ( 1/4.*( 3.*x(2)- 4.*x(4) + x(5)) .^2)))) ...

./ ( (0.3 ./( epsilon + ((13/12.*( x(3)- 2.*x(4) + x(5)).^2)+ ( 1/4.*( 3.*x(2)- 4.*x(4) + x(5)) .^2)))) ...

+ (0.6 ./( epsilon + ((13/12.*( x(2)- 2.*x(3) + x(4)).^2)+ ( 1/4.*( x(2)- x(4) ).^2)))) ...

+ (0.1 ./( epsilon + ((13/12.*( x(1) - 2.*x(2) + x(3)).^2)+ ( 1/4.*( x(1) - 4.*x(2) + 3*x(3) ).^2)))) ...

) ...

) ...

.* ( 1/3 .* x(3) + 5/6 .*x(4) - 1/6 .*x(5)) ...

+ (...

(0.6 ./( epsilon + ((13/12.*( x(2)- 2.*x(3) + x(4)).^2)+ ( 1/4.*( x(2)- x(4) ).^2)))) ...

./ ( (0.3 ./( epsilon + ((13/12.*( x(3)- 2.*x(4) + x(5)).^2)+ ( 1/4.*( 3.*x(2)- 4.*x(4) + x(5)).^2)))) ...

+ (0.6 ./( epsilon + ((13/12.*( x(2)- 2.*x(3) + x(4)).^2)+ ( 1/4.*( x(2)- x(4) ).^2)))) ...

+ (0.1 ./( epsilon + ((13/12.*( x(1) - 2.*x(2) + x(3)).^2)+ ( 1/4.*( x(1) - 4.*x(2) + 3*x(3) ).^2)))) ...

)...

) ...

.* (- 1/6 .* x(2) + 5/6 .*x(3) + 1/3 .*x(4)) ...

+ (...

(0.1 ./( epsilon + ((13/12.*( x(1) - 2.*x(2) + x(3)).^2)+ ( 1/4.*( x(1) - 4.*x(2) + 3*x(3) ).^2)))) ...

./ ( (0.3 ./( epsilon + ((13/12.*( x(3)- 2.*x(4) + x(5)).^2)+ ( 1/4.*( 3.*x(2)- 4.*x(4) + x(5)).^2)))) ...

+ (0.6 ./( epsilon + ((13/12.*( x(2)- 2.*x(3) + x(4)).^2)+ ( 1/4.*( x(2)- x(4) ).^2)))) ...

+ (0.1 ./( epsilon + ((13/12.*( x(1) - 2.*x(2) + x(3)).^2)+ ( 1/4.*( x(1) - 4.*x(2) + 3*x(3) ).^2)))) ...

)...

) ...

.* ( 1/3 .* x(1) - 7/6 .*x(2) + 11/6 .*x(3)) ...

)...

-...

( (...

(2/3 /( epsilon+(x(3) - x(2) ).^2) ) ...

./ ( (2/3 /( epsilon+(x(3) - x(2) ).^2) ) ...

+ (1/3 /( epsilon+(x(2) - x(1) ).^2) ) ...

)...

) ...

* ( 0.5 * x(2) + 0.5 * x(3)) ...

+ (...

(1/3 /( epsilon+(x(2) - x(1) ).^2) ) ...

./ ( (2/3 /( epsilon+(x(3) - x(2) ).^2) ) ...

+ (1/3 /( epsilon+(x(2) - x(1) ).^2) ) ...

)...

) ...

* (-1/2 * x(1) + 3/2 * x(2)) ...

)...

) ./ deltaZ;

%% WENO k=3

dfdt(4:Nz+2) = -u_int *( ...

( ( ...

(0.3 ./( epsilon + ((13/12.*( x(4:Nz+2)- 2.*x(5:Nz+3) + x(6:Nz+4)).^2)+ ( 1/4.*( 3.*x(3:Nz+1)- 4.*x(5:Nz+3) + x(6:Nz+4) ) .^2))).^2) ...

./ ( (0.3 ./( epsilon + ((13/12.*( x(4:Nz+2)- 2.*x(5:Nz+3) + x(6:Nz+4)).^2)+ ( 1/4.*( 3.*x(3:Nz+1)- 4.*x(5:Nz+3) + x(6:Nz+4) ) .^2))).^2) ...

+ (0.6 ./( epsilon + ((13/12.*( x(3:Nz+1)- 2.*x(4:Nz+2) + x(5:Nz+3)).^2)+ ( 1/4.*( x(3:Nz+1)- x(5:Nz+3) ).^2))).^2) ...

+ (0.1 ./( epsilon + ((13/12.*( x(2:Nz) - 2.*x(3:Nz+1) + x(4:Nz+2)).^2)+ ( 1/4.*( x(2:Nz) - 4.*x(3:Nz+1) + 3*x(4:Nz+2) ).^2))).^2) ...

) ...

) ...

.* ( 1/3 .* x(4:Nz+2) + 5/6 .*x(5:Nz+3) - 1/6 .*x(6:Nz+4)) ...

+ (...

(0.6 ./( epsilon + ((13/12.*( x(3:Nz+1)- 2.*x(4:Nz+2) + x(5:Nz+3)).^2)+ ( 1/4.*( x(3:Nz+1)- x(5:Nz+3) ).^2))).^2) ...

./ ( (0.3 ./( epsilon + ((13/12.*( x(4:Nz+2)- 2.*x(5:Nz+3) + x(6:Nz+4)).^2)+ ( 1/4.*( 3.*x(3:Nz+1)- 4.*x(5:Nz+3) + x(6:Nz+4) ).^2))).^2) ...

+ (0.6 ./( epsilon + ((13/12.*( x(3:Nz+1)- 2.*x(4:Nz+2) + x(5:Nz+3)).^2)+ ( 1/4.*( x(3:Nz+1)- x(5:Nz+3) ).^2))).^2) ...

+ (0.1 ./( epsilon + ((13/12.*( x(2:Nz) - 2.*x(3:Nz+1) + x(4:Nz+2)).^2)+ ( 1/4.*( x(2:Nz) - 4.*x(3:Nz+1) + 3*x(4:Nz+2) ).^2))).^2) ...

)...

) ...

.* (- 1/6 .* x(3:Nz+1) + 5/6 .*x(4:Nz+2) + 1/3 .*x(5:Nz+3)) ...

+ (...

(0.1 ./( epsilon + ((13/12.*( x(2:Nz) - 2.*x(3:Nz+1) + x(4:Nz+2)).^2)+ ( 1/4.*( x(2:Nz) - 4.*x(3:Nz+1) + 3*x(4:Nz+2) ).^2))).^2) ...

./ ( (0.3 ./( epsilon + ((13/12.*( x(4:Nz+2)- 2.*x(5:Nz+3) + x(6:Nz+4)).^2)+ ( 1/4.*( 3.*x(3:Nz+1)- 4.*x(5:Nz+3) + x(6:Nz+4) ).^2))).^2) ...

+ (0.6 ./( epsilon + ((13/12.*( x(3:Nz+1)- 2.*x(4:Nz+2) + x(5:Nz+3)).^2)+ ( 1/4.*( x(3:Nz+1)- x(5:Nz+3) ).^2))).^2) ...

+ (0.1 ./( epsilon + ((13/12.*( x(2:Nz) - 2.*x(3:Nz+1) + x(4:Nz+2)).^2)+ ( 1/4.*( x(2:Nz) - 4.*x(3:Nz+1) + 3*x(4:Nz+2) ).^2))).^2) ...

)...

) ...

.* ( 1/3 .* x(2:Nz) - 7/6 .*x(3:Nz+1) + 11/6 .*x(4:Nz+2)) ...

)...

-...

( ( ...

(0.3 ./( epsilon + ((13/12.*( x(3:Nz+1)- 2.*x(4:Nz+2) + x(5:Nz+3)).^2)+ ( 1/4.*( 3.*x(2:Nz) - 4.*x(4:Nz+2) + x(5:Nz+3) ) .^2))).^2) ...

./ ( (0.3 ./( epsilon + ((13/12.*( x(3:Nz+1)- 2.*x(4:Nz+2) + x(5:Nz+3)).^2)+ ( 1/4.*( 3.*x(2:Nz) - 4.*x(4:Nz+2) + x(5:Nz+3) ) .^2))).^2) ...

+ (0.6 ./( epsilon + ((13/12.*( x(2:Nz) - 2.*x(3:Nz+1) + x(4:Nz+2)).^2)+ ( 1/4.*( x(2:Nz) - x(4:Nz+2) ).^2))).^2) ...

+ (0.1 ./( epsilon + ((13/12.*( x(1:Nz-1)- 2.*x(2:Nz) + x(3:Nz+1)).^2)+ ( 1/4.*( x(1:Nz-1)- 4.*x(2:Nz) + 3*x(3:Nz+1) ).^2))).^2) ...

) ...

) ...

.* ( 1/3 .* x(3:Nz+1) + 5/6 .*x(4:Nz+2) - 1/6 .*x(5:Nz+3)) ...

+ (...

(0.6 ./( epsilon + ((13/12.*( x(2:Nz) - 2.*x(3:Nz+1) + x(4:Nz+2)).^2)+ ( 1/4.*( x(2:Nz) - x(4:Nz+2) ).^2))).^2) ...

./ ( (0.3 ./( epsilon + ((13/12.*( x(3:Nz+1)- 2.*x(4:Nz+2) + x(5:Nz+3)).^2)+ ( 1/4.*( 3.*x(2:Nz) - 4.*x(4:Nz+2) + x(5:Nz+3) ).^2))).^2) ...

+ (0.6 ./( epsilon + ((13/12.*( x(2:Nz) - 2.*x(3:Nz+1) + x(4:Nz+2)).^2)+ ( 1/4.*( x(2:Nz) - x(4:Nz+2) ).^2))).^2) ...

+ (0.1 ./( epsilon + ((13/12.*( x(1:Nz-1)- 2.*x(2:Nz) + x(3:Nz+1)).^2)+ ( 1/4.*( x(1:Nz-1)- 4.*x(2:Nz) + 3*x(3:Nz+1) ).^2))).^2) ...

)...

) ...

.* (- 1/6 .* x(2:Nz) + 5/6 .*x(3:Nz+1) + 1/3 .*x(4:Nz+2)) ...

+ (...

(0.1 ./( epsilon + ((13/12.*( x(1:Nz-1)- 2.*x(2:Nz) + x(3:Nz+1)).^2)+ ( 1/4.*( x(1:Nz-1)- 4.*x(2:Nz) + 3*x(3:Nz+1) ).^2))).^2) ...

./ ( (0.3 ./( epsilon + ((13/12.*( x(3:Nz+1)- 2.*x(4:Nz+2) + x(5:Nz+3)).^2)+ ( 1/4.*( 3.*x(2:Nz) - 4.*x(4:Nz+2) + x(5:Nz+3) ).^2))).^2) ...

+ (0.6 ./( epsilon + ((13/12.*( x(2:Nz) - 2.*x(3:Nz+1) + x(4:Nz+2)).^2)+ ( 1/4.*( x(2:Nz) - x(4:Nz+2) ).^2))).^2) ...

+ (0.1 ./( epsilon + ((13/12.*( x(1:Nz-1)- 2.*x(2:Nz) + x(3:Nz+1)).^2)+ ( 1/4.*( x(1:Nz-1)- 4.*x(2:Nz) + 3*x(3:Nz+1) ).^2))).^2) ...

)...

) ...

.* ( 1/3 .* x(1:Nz-1) - 7/6 .*x(2:Nz) + 11/6 .*x(3:Nz+1)) ...

)...

)...

./ deltaZ;

%% WENO k=2 end

dfdt(Nz+3) = -u_int * (...

( (...

(2/3 /( epsilon+(x(Nz+4) - x(Nz+3) ).^2 ).^2 ) ...

./ ( (2/3 /( epsilon+(x(Nz+4) - x(Nz+3) ).^2 ).^2 ) ...

+ (1/3 /( epsilon+(x(Nz+3) - x(Nz+2) ).^2 ).^2 ) ...

)...

) ...

* ( 0.5 * x(Nz+3) + 0.5 * x(Nz+4)) ...

+ (...

(1/3 /( epsilon+(x(Nz+3) - x(Nz+2) ).^2 ).^2 ) ...

./ ( (2/3 /( epsilon+(x(Nz+4) - x(Nz+3) ).^2 ).^2 ) ...

+ (1/3 /( epsilon+(x(Nz+3) - x(Nz+2) ).^2 ).^2 ) ...

)...

) ...

* (-1/2 * x(Nz+2) + 3/2 * x(Nz+3)) ...

) ...

- ...

( (...

(2/3 /( epsilon+(x(Nz+3) - x(Nz+2) ).^2 ).^2 ) ...

./ ( (2/3 /( epsilon+(x(Nz+3) - x(Nz+2) ).^2 ).^2 ) ...

+ (1/3 /( epsilon+(x(Nz+2) - x(Nz+1) ).^2 ).^2 ) ...

)...

) ...

* ( 0.5 * x(Nz+2) + 0.5 * x(Nz+3)) ...

+ (...

(1/3 /( epsilon+(x(Nz+2) - x(Nz+1) ).^2 ).^2 ) ...

./ ( (2/3 /( epsilon+(x(Nz+3) - x(Nz+2) ).^2 ).^2 ) ...

+ (1/3 /( epsilon+(x(Nz+2) - x(Nz+1) ).^2 ).^2 ) ...

)...

) ...

* (-1/2 * x(Nz+1) + 3/2 * x(Nz+2)) ...

)...

)./ deltaZ;

%% Upwind - end

dfdt(Nz+4) = -u_int * (x(Nz+4) - x(Nz+3)) ./ deltaZ;

end

I am getting following error:

Error using symengine>@(t,in2)reshape([-1.2e+1,(1.0./((in2(1,:)-1.0e+1).^2+1.0e-12).^2.*6.0)./(1.0./((in2(1,:)-in2(2,:)).^2+1.0e-12).^2.*(2.0./3.0)+1. .......

Too many input arguments.

Error in DBK_WENO_jacobian>@(t,x)JacFunc(u_int,deltaZ,Nz,epsilon,cFeed) (line 21)

options = odeset( 'Jacobian' , @(t,x)JacFunc(u_int, deltaZ, Nz, epsilon,cFeed), 'RelTol',1e-8, 'AbsTol',1e-8);

Error in ode15s (line 346)

dfdy = Jac(t,y,Jargs{:});

Error in DBK_WENO_jacobian (line 23)

[t,V] = ode15s(@(t,x)odeTest(t,x, u_int, deltaZ, Nz,epsilon, cFeed), tspan, v0, options);

[Torsten]

I didn't use this options command in my answer:

options = odeset( 'Jacobian' , @(t,x)JacFunc(u_int, deltaZ, Nz, epsilon,cFeed), 'RelTol',1e-8, 'AbsTol',1e-8);

[Berry]

My mistake! Thank you very much! Everything works fine now.