How to plot streamline, streakline and pathlines without using 'streamline streamline, odexx or similar functions in Matlab'

Can someone help me correct my code for tracking centre of vortex and for plotting streamlines, streaklines and pathlines by creating the code and not using Matlab 'streamline' function.

close all;
clear;
clc;

%3(a) 

time_steps = 0.01;

t=0:time_steps:100;
total_time_points=length(t);

load('particle_positions.mat')
figure
plot(xp,yp,'k.')
xlabel('X');
ylabel('Y');

for i=1:length(xp)
   
    [x_val(i),y_val(i)]=Euler_scheme(xp(i),yp(i),total_time_points,time_steps);
end

 figure
 plot(x_val,y_val,'k.')
 xlabel('X');
 ylabel('Y');
title('Particle positions')
hold on;
% Computing the center of vortex

for i=1:length(xp)/10 -1
x_center(i) =  sum(x_val(i*10:i*10+10))/10;
y_center(i) =  sum(y_val(i*10:i*10+10))/10;

end
plot(x_center,y_center,'r*')
 xlabel('X');
 ylabel('Y');


figure
plot(x_center,y_center,'r*')
 xlabel('X');
 ylabel('Y');
title('Center of Vortex');


function [x_last,y_last]=Euler_scheme(x,y,l_t,h)
 
      for j=1:l_t-1
         
x(j+1)=x(j)+(h*f(x(j),y(j)));
y(j+1)=y(j)+(h*g(x(j),y(j)));
      end
x_last=x(end);
y_last=y(end);
   end

function dxdt=f(x,y)
gamma_val = [-2 2 -2 2];
delta = 0.5;
x0 = [-1 1 -1 1];
y0 = [10 10 -10 -10];
dxdt = 0;
    for ii=1:size(x0,2)
      dxdt = dxdt - gamma_val(ii)/(2*pi)*((y-y0(1,ii))/((x-x0(1,ii))^2+(y-y0(1,ii))^2+delta^2));
    end
end
 
 function dydt=g(x,y)
gamma_val = [-2 2 -2 2];
delta = 0.5;
x0 = [-1 1 -1 1];
y0 = [10 10 -10 -10];

dydt = 0;
    for ii=1:size(x0,2)
 dydt = dydt + gamma_val(ii)/(2*pi)*((x-x0(1,ii))/((x-x0(1,ii))^2+(y-y0(1,ii))^2+delta^2));
    end
 end

Answer 1

[darova]

What about ode45? Streamline is just a trajectory, if you have vector field you can find a solution

[x,y,z] = peaks(20);
[u,v] = gradient(z);
fx = scatteredInterpolant(x(:),y(:),u(:));      % function of X velocity
fy = scatteredInterpolant(x(:),y(:),v(:));      % function of Y velocity
F = @(t,u) [fx(u(1),u(2)); fy(u(1),u(2))];      % ode45 function
quiver(x,y,u,v)
t = 0:.2:2*pi;
[x0,y0] = pol2cart(t,0.5);
for i = 1:numel(x0)
    [t,u1] = ode45(F,[0 1],[x0(i)-1.5 y0(i)+0.3]); % find solution for each initial position
    line(u1(:,1),u1(:,2),'color','r')
end