I am trying to plot a graph using mathematical formulas.

Hello everyone, I am trying to obtain the following graph using the explanations below, but I am getting a different result. I would really appreciate it if you could help me find the error.

explanations:

my code:

clear all
clc

% Parameters
epsilon0 = 8.85e-12;    
epsilon_m = 79;         
CM = 0.5;               
k_B = 1.38e-23;         
R = 1e-6;              
gamma = 1.88e-8;       
q = 1e-14;             
dt = 1e-3;             
T = 300;                
x0 = 10e-6;            
num_steps = 1000;       
N = 100000;            

epsilon = 1e-9;        
k = 1 / (4 * pi * epsilon_m);

% **Creating matrices to store FDEP force and positions**
x_dep = zeros(num_steps, N);  
x_diff = zeros(num_steps, N);
FDEP = zeros(num_steps, N);  

% **Initial Positions**
x_dep(1, :) = x0;  
x_diff(1, :) = x0;  

% **Generating random numbers for Brownian motion**
rng(0);

% **Loop for 100,000 Simulations (Each Simulation Runs Independently!)**
for sim = 1:N
    % **Generate random Gaussian noise for each particle**
    W = randn(num_steps, 1); 
    
    for i = 1:num_steps-1
        % **Calculate FDEP separately for each particle**
        FDEP(i, sim) = 4 * pi * R^3 * epsilon0 * epsilon_m * CM * (-2 * k^2 * q^2) / ((abs(x_dep(i, sim) - x0) + epsilon)^5);
        
        % **Motion under the effect of FDEP (Used for Correct Simulation)**
        x_dep(i+1, sim) = x_dep(i, sim) + (FDEP(i, sim) / gamma) * dt + sqrt(2 * k_B * T / gamma) * sqrt(dt) * W(i);
        
        % **Pure diffusion (without DEP force)**
        x_diff(i+1, sim) = x_diff(i, sim) + sqrt(2 * k_B * T / gamma) * sqrt(dt) * W(i);
    end
end

% **Calculate Mean Positions**
x_mean_dep = mean(x_dep, 2);
x_mean_diff = mean(x_diff, 2);

% **Plot the Graph**
figure;
plot((0:num_steps-1) * dt, x_mean_dep * 1e6, 'b', 'LineWidth', 1.5);
hold on;
plot((0:num_steps-1) * dt, x_mean_diff * 1e6, 'r--', 'LineWidth', 1.5);
xlabel('Time (s)');
ylabel('Position (μm)');
title('Particle Motion Under the Effect of Dielectrophoretic Force'); 
legend('Diffusion under F_{DEP}', 'Pure diffusion');
grid on;
hold off;

[Sam Chak]

Hi @Zaref Li

Didn’t you ask a similar question a few months ago? If it is indeed similar, would you like to continue the discussion there, or do you wish @Aquatris' helpful comments to be migrated here for a fresh discussion?

https://www.mathworks.com/matlabcentral/answers/2148339-how-do-i-plot-a-graph