Navier-Stokes Solver With Finite Difference Method – Unexpected Results
Dear All,
I have been working on a Matlab script to implement a Navier-Stokes solver for the lid-driven wall problem but I think there is something wrong in my math because I am not getting the expected results if I compare my plots with benchmarks I have collected from online literatures and also I am not able to infer turbolence in my results by increasing the Re number.
I have also added a little function that checks for the boundaries conditions.
The solver is using the finite difference method to update the velocity and pressure fields, it calculates the tentative velocity using the momentum equations, corrects the pressure using the Poisson equation, finally it updates the velocity by using the corrected pressure. The boundary conditions are also applied.
This is the code if any of you fancy to help!
Alex
% Navier-Stokes Solver for Incompressible, Isothermal, and Newtonian Flow
clc
% Parameters
Lx = 1; % Length in x-direction
Ly = 1; % Length in y-direction
Nx = 32; % Number of grid points in x-direction
Ny = 32; % Number of grid points in y-direction
Re = 100; % Reynolds number
U_wall = 1; % Lid-driven wall velocity
tolerance = 1e-4; % Convergence tolerance
omega = 10; % frequency of oscillations
Amp = 0.1; % amplitude of oscillations
% Discretization
dx = Lx / (Nx – 1);
dy = Ly / (Ny – 1);
x = linspace(0, Lx, Nx);
y = linspace(0, Ly, Ny);
[X, Y] = meshgrid(x, y);
% Initialize velocity and pressure fields
u = zeros(Ny, Nx);
v = zeros(Ny, Nx);
p = zeros(Ny, Nx);
% Time stepping parameters
dt = 0.001;
max_iter = 50000;
% Main loop
for iter = 1:max_iter
% Calculate lid-driven wall velocity with sinusoidal variations
U_wall = 1 + Amp*sin(omega*dt*iter);
[u, v, p] = navier_stokes_solver(u, v, p, dx, dy, dt, Re, U_wall, tolerance);
if mod(iter, 100) == 0
fprintf(‘Iteration: %dn’, iter);
pause(0.01);
end
end
%Plot Pressure Field
contourf(X, Y, p);
title(‘Pressure’);
colorbar;
pause(0.01);
% Plot velocity field
figure;
quiver(X, Y, u, v);
title(‘Velocity Field’);
xlabel(‘x’);
ylabel(‘y’);
figure;
contourf(X,Y,u’,23,’LineColor’,’none’);
title(‘U-Velocity’);
colorbar;
%Plot stream function
N = 32; xstart = max(x)*rand(N,1); ystart = max(y)*rand(N,1);
%[X,Y] = meshgrid(x,y);
figure; h=streamline(X,Y,u’,v’,xstart,ystart,[0.1, 200]);
title(‘Stream Function’); xlabel(‘x-location’); ylabel(‘y-location’)
axis(‘equal’,[0 Lx 0 Ly]); set(h,’color’,’k’)
boundary_condition_tests(u, v, U_wall);
function [u_new, v_new, p_new] = navier_stokes_solver(u, v, p, dx, dy, dt, Re, U_wall, tolerance)
[Ny, Nx] = size(u);
% Initialize temporary fields
u_new = u;
v_new = v;
p_new = p;
% Calculate the tentative velocity (u*) using momentum equations
for i = 2:Ny – 1
for j = 2:Nx – 1
% Viscous acceleration term (in x-direction)
viscous_x_term = (1 / Re) * (u(i, j + 1) – 2 * u(i, j) + u(i, j – 1)) / (dx^2);
% Viscous acceleration term (in y-direction)
viscous_y_term = (1 / Re) * (u(i + 1, j) – 2 * u(i, j) + u(i – 1, j)) / (dy^2);
u_new(i, j) = u(i, j) + dt * (viscous_x_term + viscous_y_term);
% Viscous term in x-direction
viscous_v_x_term = (1 / Re) * (v(i, j + 1) – 2 * v(i, j) + v(i, j – 1)) / (dx^2);
viscous_v_y_term = (1 / Re) * (v(i + 1, j) – 2 * v(i, j) + v(i – 1, j)) / (dy^2);
v_new(i, j) = v(i, j) + dt * (viscous_v_x_term + viscous_v_y_term);
end
end
% Pressure correction (Poisson equation)
p_initial = p_new; % Store the initial pressure field
maxe = 1;
while maxe > tolerance
p_old = p_new;
for i = 2:Ny – 1
for j = 2:Nx – 1
p_new(i, j) = (1 / 4) * (p_old(i, j + 1) + p_old(i, j – 1) …
+ p_old(i + 1, j) + p_old(i – 1, j) …
– dx * (u_new(i, j + 1) – u_new(i, j – 1)) / 2 …
– dy * (v_new(i + 1, j) – v_new(i – 1, j)) / 2);
end
end
maxe = max(max(abs(p_new – p_old)));
end
%Update velocity with pressure correction using the initial pressure field
for i = 2:Ny – 1
for j = 2:Nx – 1
u_new(i, j) = u_new(i, j) – dt * (p_initial(i, j + 1) – p_initial(i, j – 1)) / (2 * dx);
v_new(i, j) = v_new(i, j) – dt * (p_initial(i + 1, j) – p_initial(i – 1, j)) / (2 * dy);
end
end
% Boundary conditions
u_new(1, π = 0;
u_new(end, π = 0;
u_new(:, 1) = 0;
u_new(:, end) = 2 * U_wall – u_new(:, end – 1); % Moving lid
v_new(1, π = -v_new(2, :);
v_new(end, π = -v_new(end – 1, :);
v_new(:, 1) = 0;
v_new(:, end) = 0;
p_new(1, π = p_new(2, :);
p_new(end, π = p_new(end – 1, :);
p_new(:, 1) = p_new(:, 2);
p_new(:, end) = p_new(:, end – 1);
end
function boundary_condition_tests(u, v, U_wall)
fprintf(‘Boundary Condition Testn’);
% Bottom wall u-velocity
bottom_wall_u = u(1, :);
bottom_wall_u_test = all(abs(bottom_wall_u) < 1e-6);
if bottom_wall_u_test
fprintf(‘Bottom wall u-velocity: Passn’);
else
fprintf(‘Bottom wall u-velocity: Fail (Actual value: %s)n’, mat2str(bottom_wall_u));
end
% Top wall u-velocity
top_wall_u = u(end, :);
top_wall_u_test = all(abs(top_wall_u – 2 * U_wall) < 1e-6);
if top_wall_u_test
fprintf(‘Top wall u-velocity: Passn’);
else
fprintf(‘Top wall u-velocity: Fail (Actual value: %s)n’, mat2str(top_wall_u));
end
% Left and right wall u-velocity
left_wall_u = u(:, 1);
right_wall_u = u(:, end);
left_right_wall_u_test = all(abs(left_wall_u) < 1e-6) && all(abs(right_wall_u) < 1e-6);
if left_right_wall_u_test
fprintf(‘Left and right wall u-velocity: Passn’);
else
fprintf(‘Left and right wall u-velocity: Fail (Actual value: Left – %s, Right – %s)n’, mat2str(left_wall_u), mat2str(right_wall_u));
end
% Left and right wall v-velocity
left_wall_v = v(:, 1);
right_wall_v = v(:, end);
left_right_wall_v_test = all(abs(left_wall_v) < 1e-6) && all(abs(right_wall_v) < 1e-6);
if left_right_wall_v_test
fprintf(‘Left and right wall v-velocity: Passn’);
else
fprintf(‘Left and right wall v-velocity: Fail (Actual value: Left – %s, Right – %s)n’, mat2str(left_wall_v), mat2str(right_wall_v));
end
% Bottom and top wall v-velocity
bottom_wall_v = v(1, :);
top_wall_v = v(end, :);
bottom_top_wall_v_test = all(abs(bottom_wall_v + v(2, :)) < 1e-6) && all(abs(top_wall_v + v(end – 1, :)) < 1e-6);
if bottom_top_wall_v_test
fprintf(‘Bottom and top wall v-velocity: Passn’);
else
fprintf(‘Bottom and top wall v-velocity: Fail (Actual value: Bottom – %s, Top – %s)n’, mat2str(bottom_wall_v), mat2str(top_wall_v));
end
% Count the total failed boundary conditions
failed_count = sum([~bottom_wall_u_test, ~top_wall_u_test, ~left_right_wall_u_test, ~left_right_wall_v_test, ~bottom_top_wall_v_test]);
fprintf(‘Total failed boundary conditions: %dn’, failed_count);
endDear All,
I have been working on a Matlab script to implement a Navier-Stokes solver for the lid-driven wall problem but I think there is something wrong in my math because I am not getting the expected results if I compare my plots with benchmarks I have collected from online literatures and also I am not able to infer turbolence in my results by increasing the Re number.
I have also added a little function that checks for the boundaries conditions.
The solver is using the finite difference method to update the velocity and pressure fields, it calculates the tentative velocity using the momentum equations, corrects the pressure using the Poisson equation, finally it updates the velocity by using the corrected pressure. The boundary conditions are also applied.
This is the code if any of you fancy to help!
Alex
% Navier-Stokes Solver for Incompressible, Isothermal, and Newtonian Flow
clc
% Parameters
Lx = 1; % Length in x-direction
Ly = 1; % Length in y-direction
Nx = 32; % Number of grid points in x-direction
Ny = 32; % Number of grid points in y-direction
Re = 100; % Reynolds number
U_wall = 1; % Lid-driven wall velocity
tolerance = 1e-4; % Convergence tolerance
omega = 10; % frequency of oscillations
Amp = 0.1; % amplitude of oscillations
% Discretization
dx = Lx / (Nx – 1);
dy = Ly / (Ny – 1);
x = linspace(0, Lx, Nx);
y = linspace(0, Ly, Ny);
[X, Y] = meshgrid(x, y);
% Initialize velocity and pressure fields
u = zeros(Ny, Nx);
v = zeros(Ny, Nx);
p = zeros(Ny, Nx);
% Time stepping parameters
dt = 0.001;
max_iter = 50000;
% Main loop
for iter = 1:max_iter
% Calculate lid-driven wall velocity with sinusoidal variations
U_wall = 1 + Amp*sin(omega*dt*iter);
[u, v, p] = navier_stokes_solver(u, v, p, dx, dy, dt, Re, U_wall, tolerance);
if mod(iter, 100) == 0
fprintf(‘Iteration: %dn’, iter);
pause(0.01);
end
end
%Plot Pressure Field
contourf(X, Y, p);
title(‘Pressure’);
colorbar;
pause(0.01);
% Plot velocity field
figure;
quiver(X, Y, u, v);
title(‘Velocity Field’);
xlabel(‘x’);
ylabel(‘y’);
figure;
contourf(X,Y,u’,23,’LineColor’,’none’);
title(‘U-Velocity’);
colorbar;
%Plot stream function
N = 32; xstart = max(x)*rand(N,1); ystart = max(y)*rand(N,1);
%[X,Y] = meshgrid(x,y);
figure; h=streamline(X,Y,u’,v’,xstart,ystart,[0.1, 200]);
title(‘Stream Function’); xlabel(‘x-location’); ylabel(‘y-location’)
axis(‘equal’,[0 Lx 0 Ly]); set(h,’color’,’k’)
boundary_condition_tests(u, v, U_wall);
function [u_new, v_new, p_new] = navier_stokes_solver(u, v, p, dx, dy, dt, Re, U_wall, tolerance)
[Ny, Nx] = size(u);
% Initialize temporary fields
u_new = u;
v_new = v;
p_new = p;
% Calculate the tentative velocity (u*) using momentum equations
for i = 2:Ny – 1
for j = 2:Nx – 1
% Viscous acceleration term (in x-direction)
viscous_x_term = (1 / Re) * (u(i, j + 1) – 2 * u(i, j) + u(i, j – 1)) / (dx^2);
% Viscous acceleration term (in y-direction)
viscous_y_term = (1 / Re) * (u(i + 1, j) – 2 * u(i, j) + u(i – 1, j)) / (dy^2);
u_new(i, j) = u(i, j) + dt * (viscous_x_term + viscous_y_term);
% Viscous term in x-direction
viscous_v_x_term = (1 / Re) * (v(i, j + 1) – 2 * v(i, j) + v(i, j – 1)) / (dx^2);
viscous_v_y_term = (1 / Re) * (v(i + 1, j) – 2 * v(i, j) + v(i – 1, j)) / (dy^2);
v_new(i, j) = v(i, j) + dt * (viscous_v_x_term + viscous_v_y_term);
end
end
% Pressure correction (Poisson equation)
p_initial = p_new; % Store the initial pressure field
maxe = 1;
while maxe > tolerance
p_old = p_new;
for i = 2:Ny – 1
for j = 2:Nx – 1
p_new(i, j) = (1 / 4) * (p_old(i, j + 1) + p_old(i, j – 1) …
+ p_old(i + 1, j) + p_old(i – 1, j) …
– dx * (u_new(i, j + 1) – u_new(i, j – 1)) / 2 …
– dy * (v_new(i + 1, j) – v_new(i – 1, j)) / 2);
end
end
maxe = max(max(abs(p_new – p_old)));
end
%Update velocity with pressure correction using the initial pressure field
for i = 2:Ny – 1
for j = 2:Nx – 1
u_new(i, j) = u_new(i, j) – dt * (p_initial(i, j + 1) – p_initial(i, j – 1)) / (2 * dx);
v_new(i, j) = v_new(i, j) – dt * (p_initial(i + 1, j) – p_initial(i – 1, j)) / (2 * dy);
end
end
% Boundary conditions
u_new(1, π = 0;
u_new(end, π = 0;
u_new(:, 1) = 0;
u_new(:, end) = 2 * U_wall – u_new(:, end – 1); % Moving lid
v_new(1, π = -v_new(2, :);
v_new(end, π = -v_new(end – 1, :);
v_new(:, 1) = 0;
v_new(:, end) = 0;
p_new(1, π = p_new(2, :);
p_new(end, π = p_new(end – 1, :);
p_new(:, 1) = p_new(:, 2);
p_new(:, end) = p_new(:, end – 1);
end
function boundary_condition_tests(u, v, U_wall)
fprintf(‘Boundary Condition Testn’);
% Bottom wall u-velocity
bottom_wall_u = u(1, :);
bottom_wall_u_test = all(abs(bottom_wall_u) < 1e-6);
if bottom_wall_u_test
fprintf(‘Bottom wall u-velocity: Passn’);
else
fprintf(‘Bottom wall u-velocity: Fail (Actual value: %s)n’, mat2str(bottom_wall_u));
end
% Top wall u-velocity
top_wall_u = u(end, :);
top_wall_u_test = all(abs(top_wall_u – 2 * U_wall) < 1e-6);
if top_wall_u_test
fprintf(‘Top wall u-velocity: Passn’);
else
fprintf(‘Top wall u-velocity: Fail (Actual value: %s)n’, mat2str(top_wall_u));
end
% Left and right wall u-velocity
left_wall_u = u(:, 1);
right_wall_u = u(:, end);
left_right_wall_u_test = all(abs(left_wall_u) < 1e-6) && all(abs(right_wall_u) < 1e-6);
if left_right_wall_u_test
fprintf(‘Left and right wall u-velocity: Passn’);
else
fprintf(‘Left and right wall u-velocity: Fail (Actual value: Left – %s, Right – %s)n’, mat2str(left_wall_u), mat2str(right_wall_u));
end
% Left and right wall v-velocity
left_wall_v = v(:, 1);
right_wall_v = v(:, end);
left_right_wall_v_test = all(abs(left_wall_v) < 1e-6) && all(abs(right_wall_v) < 1e-6);
if left_right_wall_v_test
fprintf(‘Left and right wall v-velocity: Passn’);
else
fprintf(‘Left and right wall v-velocity: Fail (Actual value: Left – %s, Right – %s)n’, mat2str(left_wall_v), mat2str(right_wall_v));
end
% Bottom and top wall v-velocity
bottom_wall_v = v(1, :);
top_wall_v = v(end, :);
bottom_top_wall_v_test = all(abs(bottom_wall_v + v(2, :)) < 1e-6) && all(abs(top_wall_v + v(end – 1, :)) < 1e-6);
if bottom_top_wall_v_test
fprintf(‘Bottom and top wall v-velocity: Passn’);
else
fprintf(‘Bottom and top wall v-velocity: Fail (Actual value: Bottom – %s, Top – %s)n’, mat2str(bottom_wall_v), mat2str(top_wall_v));
end
% Count the total failed boundary conditions
failed_count = sum([~bottom_wall_u_test, ~top_wall_u_test, ~left_right_wall_u_test, ~left_right_wall_v_test, ~bottom_top_wall_v_test]);
fprintf(‘Total failed boundary conditions: %dn’, failed_count);
endΒ Dear All,
I have been working on a Matlab script to implement a Navier-Stokes solver for the lid-driven wall problem but I think there is something wrong in my math because I am not getting the expected results if I compare my plots with benchmarks I have collected from online literatures and also I am not able to infer turbolence in my results by increasing the Re number.
I have also added a little function that checks for the boundaries conditions.
The solver is using the finite difference method to update the velocity and pressure fields, it calculates the tentative velocity using the momentum equations, corrects the pressure using the Poisson equation, finally it updates the velocity by using the corrected pressure. The boundary conditions are also applied.
This is the code if any of you fancy to help!
Alex
% Navier-Stokes Solver for Incompressible, Isothermal, and Newtonian Flow
clc
% Parameters
Lx = 1; % Length in x-direction
Ly = 1; % Length in y-direction
Nx = 32; % Number of grid points in x-direction
Ny = 32; % Number of grid points in y-direction
Re = 100; % Reynolds number
U_wall = 1; % Lid-driven wall velocity
tolerance = 1e-4; % Convergence tolerance
omega = 10; % frequency of oscillations
Amp = 0.1; % amplitude of oscillations
% Discretization
dx = Lx / (Nx – 1);
dy = Ly / (Ny – 1);
x = linspace(0, Lx, Nx);
y = linspace(0, Ly, Ny);
[X, Y] = meshgrid(x, y);
% Initialize velocity and pressure fields
u = zeros(Ny, Nx);
v = zeros(Ny, Nx);
p = zeros(Ny, Nx);
% Time stepping parameters
dt = 0.001;
max_iter = 50000;
% Main loop
for iter = 1:max_iter
% Calculate lid-driven wall velocity with sinusoidal variations
U_wall = 1 + Amp*sin(omega*dt*iter);
[u, v, p] = navier_stokes_solver(u, v, p, dx, dy, dt, Re, U_wall, tolerance);
if mod(iter, 100) == 0
fprintf(‘Iteration: %dn’, iter);
pause(0.01);
end
end
%Plot Pressure Field
contourf(X, Y, p);
title(‘Pressure’);
colorbar;
pause(0.01);
% Plot velocity field
figure;
quiver(X, Y, u, v);
title(‘Velocity Field’);
xlabel(‘x’);
ylabel(‘y’);
figure;
contourf(X,Y,u’,23,’LineColor’,’none’);
title(‘U-Velocity’);
colorbar;
%Plot stream function
N = 32; xstart = max(x)*rand(N,1); ystart = max(y)*rand(N,1);
%[X,Y] = meshgrid(x,y);
figure; h=streamline(X,Y,u’,v’,xstart,ystart,[0.1, 200]);
title(‘Stream Function’); xlabel(‘x-location’); ylabel(‘y-location’)
axis(‘equal’,[0 Lx 0 Ly]); set(h,’color’,’k’)
boundary_condition_tests(u, v, U_wall);
function [u_new, v_new, p_new] = navier_stokes_solver(u, v, p, dx, dy, dt, Re, U_wall, tolerance)
[Ny, Nx] = size(u);
% Initialize temporary fields
u_new = u;
v_new = v;
p_new = p;
% Calculate the tentative velocity (u*) using momentum equations
for i = 2:Ny – 1
for j = 2:Nx – 1
% Viscous acceleration term (in x-direction)
viscous_x_term = (1 / Re) * (u(i, j + 1) – 2 * u(i, j) + u(i, j – 1)) / (dx^2);
% Viscous acceleration term (in y-direction)
viscous_y_term = (1 / Re) * (u(i + 1, j) – 2 * u(i, j) + u(i – 1, j)) / (dy^2);
u_new(i, j) = u(i, j) + dt * (viscous_x_term + viscous_y_term);
% Viscous term in x-direction
viscous_v_x_term = (1 / Re) * (v(i, j + 1) – 2 * v(i, j) + v(i, j – 1)) / (dx^2);
viscous_v_y_term = (1 / Re) * (v(i + 1, j) – 2 * v(i, j) + v(i – 1, j)) / (dy^2);
v_new(i, j) = v(i, j) + dt * (viscous_v_x_term + viscous_v_y_term);
end
end
% Pressure correction (Poisson equation)
p_initial = p_new; % Store the initial pressure field
maxe = 1;
while maxe > tolerance
p_old = p_new;
for i = 2:Ny – 1
for j = 2:Nx – 1
p_new(i, j) = (1 / 4) * (p_old(i, j + 1) + p_old(i, j – 1) …
+ p_old(i + 1, j) + p_old(i – 1, j) …
– dx * (u_new(i, j + 1) – u_new(i, j – 1)) / 2 …
– dy * (v_new(i + 1, j) – v_new(i – 1, j)) / 2);
end
end
maxe = max(max(abs(p_new – p_old)));
end
%Update velocity with pressure correction using the initial pressure field
for i = 2:Ny – 1
for j = 2:Nx – 1
u_new(i, j) = u_new(i, j) – dt * (p_initial(i, j + 1) – p_initial(i, j – 1)) / (2 * dx);
v_new(i, j) = v_new(i, j) – dt * (p_initial(i + 1, j) – p_initial(i – 1, j)) / (2 * dy);
end
end
% Boundary conditions
u_new(1, π = 0;
u_new(end, π = 0;
u_new(:, 1) = 0;
u_new(:, end) = 2 * U_wall – u_new(:, end – 1); % Moving lid
v_new(1, π = -v_new(2, :);
v_new(end, π = -v_new(end – 1, :);
v_new(:, 1) = 0;
v_new(:, end) = 0;
p_new(1, π = p_new(2, :);
p_new(end, π = p_new(end – 1, :);
p_new(:, 1) = p_new(:, 2);
p_new(:, end) = p_new(:, end – 1);
end
function boundary_condition_tests(u, v, U_wall)
fprintf(‘Boundary Condition Testn’);
% Bottom wall u-velocity
bottom_wall_u = u(1, :);
bottom_wall_u_test = all(abs(bottom_wall_u) < 1e-6);
if bottom_wall_u_test
fprintf(‘Bottom wall u-velocity: Passn’);
else
fprintf(‘Bottom wall u-velocity: Fail (Actual value: %s)n’, mat2str(bottom_wall_u));
end
% Top wall u-velocity
top_wall_u = u(end, :);
top_wall_u_test = all(abs(top_wall_u – 2 * U_wall) < 1e-6);
if top_wall_u_test
fprintf(‘Top wall u-velocity: Passn’);
else
fprintf(‘Top wall u-velocity: Fail (Actual value: %s)n’, mat2str(top_wall_u));
end
% Left and right wall u-velocity
left_wall_u = u(:, 1);
right_wall_u = u(:, end);
left_right_wall_u_test = all(abs(left_wall_u) < 1e-6) && all(abs(right_wall_u) < 1e-6);
if left_right_wall_u_test
fprintf(‘Left and right wall u-velocity: Passn’);
else
fprintf(‘Left and right wall u-velocity: Fail (Actual value: Left – %s, Right – %s)n’, mat2str(left_wall_u), mat2str(right_wall_u));
end
% Left and right wall v-velocity
left_wall_v = v(:, 1);
right_wall_v = v(:, end);
left_right_wall_v_test = all(abs(left_wall_v) < 1e-6) && all(abs(right_wall_v) < 1e-6);
if left_right_wall_v_test
fprintf(‘Left and right wall v-velocity: Passn’);
else
fprintf(‘Left and right wall v-velocity: Fail (Actual value: Left – %s, Right – %s)n’, mat2str(left_wall_v), mat2str(right_wall_v));
end
% Bottom and top wall v-velocity
bottom_wall_v = v(1, :);
top_wall_v = v(end, :);
bottom_top_wall_v_test = all(abs(bottom_wall_v + v(2, :)) < 1e-6) && all(abs(top_wall_v + v(end – 1, :)) < 1e-6);
if bottom_top_wall_v_test
fprintf(‘Bottom and top wall v-velocity: Passn’);
else
fprintf(‘Bottom and top wall v-velocity: Fail (Actual value: Bottom – %s, Top – %s)n’, mat2str(bottom_wall_v), mat2str(top_wall_v));
end
% Count the total failed boundary conditions
failed_count = sum([~bottom_wall_u_test, ~top_wall_u_test, ~left_right_wall_u_test, ~left_right_wall_v_test, ~bottom_top_wall_v_test]);
fprintf(‘Total failed boundary conditions: %dn’, failed_count);
endΒ navier stokes solver, finite differencesΒ MATLAB Answers β New Questions
β