Numerical differential equation solver for initial conditions not at y0

I've been trying to convert a part of mathematica code to a Matlab one, in order to integrate it for other calculations, but I am having some trouble getting the solver to run.

If by any chance this helps, this is the script that I am trying to convert:

And this is my approach.

syms Uf(t) delta y

d = 0.2000;
k_Karman = 0.4;
sol_Um = 3.2947;
omega = 0.5405;
Tp = 29.0823;

ks = 2.5 * d *10^-3  % [m]
k_Karman = 0.4      % von Karman constant

U_R04 = sol_Um * sin(omega*t)
U0_R04 = Uf/k_Karman*log(30*y/ks)

eq_R04 = sol_Um * sin(omega*t) == Uf/k_Karman*log(30*delta/ks)
sol_delta = solve(eq_R04,delta) % solve for delta

ME_R04 = int(diff(U0_R04,t),y,[0 sol_delta]) == sol_delta*diff(U_R04,t)-Uf^2

I believe that this far, I would be getting a matching result from the example below, however, I am unable to prepare the command line for the evaluation.

This was my initial attempt, but I have a feeling that is the wrong approach to take for this specific case as I am getting an error for the 'Vars' definition.

[VF,Sbs] = odeToVectorField(ME_R04);
Uffcn = matlabFunction(VF, 'Vars',{x,Y});
[x,Ufx] = ode15s(Uffcn, [0 Tp/2], 0);

Error using sym/matlabFunction>getOptions

The value of 'Vars' is invalid. 'Vars' value must be a character vector, a 1-dimensional cell array of character vectors, a 1-dimensional cell array of symbolic variables

or arrays of symbolic variables, or an array of symbolic variables.

Initial attempt was done for initial conditions of y0=0, as I am not aware of a way to define it the way mathematica script allows for Uf(0.001*Tp)=0.001*sol_Um

[Rihards]

Edited the last part as follows:

[VF,Sbs] = odeToVectorField(ME_R04);
Uffcn = matlabFunction(VF, 'Vars',{t,Y});
init_cond = 0.0032947071555424111234342753369588 % Gives error when defined as 0.001*Tp
[t,Ufx] = ode15s(Uffcn, [0.001*Tp Tp/2],init_cond);

figure(9);
plot(t, Ufx,'b','linewidth',1)

And I am curious if anyone can explain, why it solves here, but gives an error when the init_cond is defined as

init_cond = 0.001*Tp

Answer 1

[Torsten]

The solution seems to be very sensitive to the initial condition.

syms Uf(t) delta y
d = 0.2000;
k_Karman = 0.4;
sol_Um = 3.2947;
omega = 0.5405;
Tp = 29.0823;
ks = 2.5 * d *10^-3;  % [m]
k_Karman = 0.4;      % von Karman constant
U_R04 = sol_Um * sin(omega*t);
U0_R04 = Uf/k_Karman*log(30*y/ks);
eq_R04 = sol_Um * sin(omega*t) == Uf/k_Karman*log(30*delta/ks);
sol_delta = solve(eq_R04,delta); % solve for delta

Warning: Solutions are only valid under certain conditions. To include parameters and conditions in the solution, specify the 'ReturnConditions' value as 'true'.

ME_R04 = int(diff(U0_R04,t),y,[0 sol_delta]) == sol_delta*diff(U_R04,t)-Uf^2;
[VF,Sbs] = odeToVectorField(ME_R04);
Uffcn = matlabFunction(VF);
init_cond = 0.02;%0.001*Tp%0.0032947071555424111234342753369588; % Gives error when defined as 0.001*Tp
options = odeset('RelTol',1e-8,'AbsTol',1e-8);
[t,Ufx] = ode15s(@(t,y)Uffcn(y,t), [0.001*Tp 5.8],init_cond,options);
figure(1);
plot(t, Ufx,'b','linewidth',1)

[Rihards]

I am getting very similar result with this. Is there a way to increase the number of solutions for the plot? Currently I got 72 values for t and Ufx, but I believe the very peak value is somewhere in between them.

syms Uf(t) delta y

ks = 2.5 * d *10^-3  % [m]
k_Karman = 0.4      % von Karman constant

U_R04 = sol_Um * sin(omega*t)
U0_R04 = Uf/k_Karman*log(30*y/ks)

eq_R04 = sol_Um * sin(omega*t) == Uf/k_Karman*log(30*delta/ks)
sol_delta = solve(eq_R04,delta) % solve for delta

ME_R04 = int(diff(U0_R04,t),y,[0 sol_delta]) == sol_delta*diff(U_R04,t)-Uf^2

[VF,Sbs] = odeToVectorField(ME_R04);
Uffcn = matlabFunction(VF, 'Vars',{t,Y});
init_cond = 0.0032947071555424111234342753369588 % Gives error when defined directly as 0.001*sol_Um
[t,Ufx] = ode15s(Uffcn, [0.001*T T/2],init_cond);

figure(9);
plot(t, Ufx,'b','linewidth',1)
title('R04-Task2','interpreter','latex','fontsize',16)
xlabel('Time, [s]','interpreter','latex','fontsize',16)
ylabel('Friction velocity Uf, [m/s]','interpreter','latex','fontsize',16)
grid on

[Torsten]

nt = 150;
[t,Ufx] = ode15s(Uffcn, linspace(0.001*T T/2,nt),init_cond);

instead of

[t,Ufx] = ode15s(Uffcn, [0.001*T T/2],init_cond);