How to get the optimization of the object function with steepest descent method?

suppose that：

 %
E1=exp(-(x1-mul1)^2/(2*sigma1^2));
E2=exp(-(x2-mul2)^2/(2*sigma2^2));
……
En=exp(-(xn-muln)^2/(2*sigman^2));

Object function :f=E1*E2*…*En;

input={x1,x2,…,xn}=[1,2,3,4,5];

Initial value of parameter mul and sigma as following:

ml={mul1,mul2,…,muln}=[2,3,4,5,6];

sa={sigma1,sigma2,…,sigman}=[1,2,3,4,5];

There is the code i have written, existing some problem during running,looking forward your help!

%  
clear all;
ml=[2,3,4,5,6];%initial value of mul
sa=[1,2,3,4,5];%initial value of sigma
input=[1,2,3,4,5];
[m,N]=size(input);
var=[];
mul=[];
sigma=[];
lambda=[]; %step size of the object function
mul=sym(mul);
sigma=sym(sigma);
lambda=sym(lambda);
for i=1:N
   mul(i)=['mul',num2str(i)];
   sigma(i)=['sigma',num2str(i)];
   lambda(i)=['lambda',num2str(i)];
end
var=transpose([mul,sigma]);
epoch=10; %the number of iteration 
mem=[];
mem=exp(-(input-mul).^2./(2*sigma.^2))
temp=prod(mem) % object function
esp=0.001;
mulgrad=[];
sigmagrad=[];
for i=1:epoch
   for j=1:N
       mulgrad(i)=diff(temp,(mul(i))) %derivative of mul
       mulgrad(i)=subs(mulgrad(i),mul(i),ml(i))
       sigmagrad(i)=diff(temp,sigma(i))
       sigmagrad(i)=subs(sigmagrad(i),sigma(i),sa(i))
   end
   norm=sum(sqrt(mulgrad^2+sigmagrad^2))/N;
   if norm<0.001
       printf('result= %g',temp);
       break;
   end 
   mul=mul+lambda.*mulgrad
   sigma=sigma+lambda.*sigmagrads
   mem=exp(-(input-mul).^2./(2*sigma.^2))
  temp=prod(mem)
  for j=1:epoch
  d(j)=diff(temp,lambda(i))
  end
%  how to get the step size ?we need update N step-size simultaneously,it
%  is different for me.
  lambdaNew=solve(d,'lambda')
  mul=mul+lambdaNew.*mul1
  sigma=sigma+lambdaNew.*sigma1
  mem=exp(-(input-mul).^2./(2*sigma.^2))
  temp=prod(mem)
end
end