syms x
L = 2*pi;
s = rand;
si = L*s;
G5_6 = gamma(5/6); G1_3 = gamma(1/3); k_di = 10000/35;
C = (55/9)*(G5_6/(sqrt(pi)*G1_3));
sigma = 2; ke = 5; l = (2*pi)/5;
n = 1; N=10; u = 0; eta = l/(10^4);
nke = n/ke; k_eta = (2*pi)/eta;
sigma_n = k_eta/10; a_n = 0;
for i = n:N
E_vk = C*((sigma^2)/ke)*((nke)^4)/((1+nke^2)^(17/6));
E_n = E_vk*exp((-(n-k_di)^2)/(sigma_n^2));
a_n = ((2*E_n)^(1/2)) + a_n;
u = (a_n*sin(n*x+si) + u);
end
xnum = linspace(0,10*pi,1000);
unum = subs(u,x,xnum);
plot(xnum,unum)
