Numerical Integration in matlab

what is reliable way to perform a integartion numerically in matlab? I would like to get numerical result after integration of a large function .

Is that Trapezoidal rule or Simpson's rule is reliable to do that? Pl somebody tell me what is the best way for that.

Otherwise, when I use 'int()' command directly, Matlab takes very very huge time almost 5-6 hours.

[James Tursa]

Please give us the details on what you are trying to integrate, and over what range.

[Star Strider]

Reference: 3D Plotting in matlab

[AVM]

@James: Thanks . Pl see this.

clc
close all 

syms 'theta' 'phi' g 
g=1;
thet=4;
h=[ cos(theta) 0 sin(theta)*exp(-1i*phi) 0
    0 cos(theta) g sin(theta)*exp(-1i*phi)
    sin(theta)*exp(1i*phi) g -cos(theta) 0
    0 sin(theta)*exp(1i*phi) 0 -cos(theta)  ];

a1=(exp(-phi*2i)*(cos(theta)^2 + sin(theta)^2 + g^2/2 - (g*(4*sin(theta)^2 + g^2)^(1/2))/2)^(3/2))/(g*sin(theta)^2) - (exp(-phi*2i)*cos(theta)*(g^2 + cos(theta)^2 + sin(theta)^2))/(g*sin(theta)^2) + (exp(-phi*2i)*cos(theta)*(cos(theta)^2 + sin(theta)^2 + g^2/2 - (g*(4*sin(theta)^2 + g^2)^(1/2))/2))/(g*sin(theta)^2) - (exp(-phi*2i)*(g^2 + cos(theta)^2 + sin(theta)^2)*(cos(theta)^2 + sin(theta)^2 + g^2/2 - (g*(4*sin(theta)^2 + g^2)^(1/2))/2)^(1/2))/(g*sin(theta)^2);
b1=(exp(-phi*1i)*(cos(theta)^2 + sin(theta)^2 + g^2/2 - (g*(4*sin(theta)^2 + g^2)^(1/2))/2)^(1/2))/sin(theta) + (exp(-phi*1i)*cos(theta))/sin(theta);
c1=-(exp(-phi*1i)*(cos(theta)^2 + sin(theta)^2))/(g*sin(theta)) + (exp(-phi*1i)*(cos(theta)^2 + sin(theta)^2 + g^2/2 - (g*(4*sin(theta)^2 + g^2)^(1/2))/2))/(g*sin(theta));
d1=1;

m=sqrt(a1*conj(a1)+b1*conj(b1)+c1*conj(c1)+d1*conj(d1));

u=1/m*[ a1
        b1
        c1
        d1 ];
    
v=diff(u,phi);
r=dot(u,v);

f(theta,g)=1i/pi*int(r,phi,0,2*pi); %% This integration taking so much time

f = simplify(f, 'Steps',500);
ffcn = matlabFunction(f);

theta = linspace(0.001,4, 30);
g = linspace(0.001,10, 30);
[Th,G] = meshgrid(theta, g);

F=ffcn(Th,G);

%max(imag(F(:)))
%min(imag(F(:)))
figure
mesh(Th, G, F)
colormap(cool)
grid on
xlabel('\bf\theta','FontSize',14)
ylabel('\bf\alpha','FontSize',14)
zlabel('\bf\itf','FontSize',14)

%% This entire process takes long time, I quit that in between.

[Walter Roberson]

simplify( r) before doing the integration; it compacts down quite a bit. Also,

assume(phi>=0)
assume(theta>=0)

before doing the simplify()

Answer 1

[Walter Roberson]

There is no reliable numeric integration method in any programming language. For any given numeric integration method, it is possible to construct a function which the numeric integration method will return an arbitrarily wrong solution.

Trapazoidal Rule and Simpson's Rule are not reliable at all.

You could consider using integral(); provide AbsTol and RelTol parameters, and provide WayPoints whenever meaningful.

You could consider using vpaintegral() with similar parameters.

[Walter Roberson]

Once you

assume(phi >= 0)

which we can do because you integrate r over phi from 0 to 2*Pi

then when you simplify(r ) the phi drops out and you end up with an expression that is only in theta. Because it is constant in the variable of integration, the integral is just the difference in integration bounds times the value of the simplified expression.

[AVM]

@ Walter: Can I take steps size smaller than 500? Bcz it is taking so much time for simplification r.

[Walter Roberson]

>> assume(phi>=0); assume(theta>=0);
>> tic;symvar(simplify(r)),toc
 
ans =
 
theta
 
Elapsed time is 0.389983 seconds.

And therefore there is no need to int() with respect to phi, since phi dropped out of the expression.

[AVM]

@Walter: Thanks a lot. The plotting is beautifully coming. I am really greatful to you.

[AVM]

@Walter: Is there any way to call the particular eigen vector of a matrix? Actually, in my current programming, this a1,b1,c1,d1 are basically components of a particular eigen vector of h matirx. I just paste them from command palte. But I would like to avoid that big expression. Is there any way to call it directly?

[Walter Roberson]

Eigenvectors are not unique so it is not clear what a particular eigenvector is.

If you use the two output form of eig() then some of the eigenvectors appear as columns of the first output. If the eigenvalues are unique then the eigenvectors output span the space.

[AVM]

@walter: Thanks ..that I know. Here I attach my code, Why this is taking so much time to reveal the eigen values of the finalrho? As if the process is never ending. You Pl look my code.

clc
close all
syms 'theta' 'alpha' 'phi' k
%% sigma matrices are here
sigma1=[0 1;1 0];
sigma2=[0 -1i;1i 0];
sigma3=[1 0;0 -1];
sigmap=1/2*(sigma1+1i*sigma2);
sigmam=1/2*(sigma1-1i*sigma2);
%% unit vetor on sphere
n=[sin(theta)*cos(phi);sin(theta)*sin(phi);cos(theta)];

identity1=eye(2);
identity2=identity1;

p1=sigma1*sin(theta)*cos(phi)+sigma2*sin(theta)*sin(phi)+sigma3*cos(theta);
p2=kron(sigmap,sigmap);
p3=kron(sigmam,sigmam);

h=1/2*alpha*kron(p1,identity2)+k*(p2+p3);

[V,L]=eig(h);

%% The componens of one of the eigen vectors of h above pasted from command palte

a1=(4*((alpha^2*cos(theta)^2)/4 - (k*(k^2 + alpha^2*cos(phi)^2*sin(theta)^2 + alpha^2*sin(phi)^2*sin(theta)^2)^(1/2))/2 + k^2/2 + (alpha^2*cos(phi)^2*sin(theta)^2)/4 + (alpha^2*sin(phi)^2*sin(theta)^2)/4)^(3/2))/(alpha^2*k*cos(phi)^2*sin(theta)^2 + alpha^2*k*sin(phi)^2*sin(theta)^2) - (4*k^2*cos(theta) + alpha^2*cos(theta)^3 + alpha^2*cos(phi)^2*cos(theta)*sin(theta)^2 + alpha^2*cos(theta)*sin(phi)^2*sin(theta)^2)/(2*(alpha*k*cos(phi)^2*sin(theta)^2 + alpha*k*sin(phi)^2*sin(theta)^2)) + (2*cos(theta)*((alpha^2*cos(theta)^2)/4 - (k*(k^2 + alpha^2*cos(phi)^2*sin(theta)^2 + alpha^2*sin(phi)^2*sin(theta)^2)^(1/2))/2 + k^2/2 + (alpha^2*cos(phi)^2*sin(theta)^2)/4 + (alpha^2*sin(phi)^2*sin(theta)^2)/4))/(alpha*k*cos(phi)^2*sin(theta)^2 + alpha*k*sin(phi)^2*sin(theta)^2) - ((alpha^2*cos(theta)^2 + 4*k^2 + alpha^2*cos(phi)^2*sin(theta)^2 + alpha^2*sin(phi)^2*sin(theta)^2)*((alpha^2*cos(theta)^2)/4 - (k*(k^2 + alpha^2*cos(phi)^2*sin(theta)^2 + alpha^2*sin(phi)^2*sin(theta)^2)^(1/2))/2 + k^2/2 + (alpha^2*cos(phi)^2*sin(theta)^2)/4 + (alpha^2*sin(phi)^2*sin(theta)^2)/4)^(1/2))/(alpha^2*k*cos(phi)^2*sin(theta)^2 + alpha^2*k*sin(phi)^2*sin(theta)^2);
b1=-(((alpha^2*cos(theta)^2)/4 - (k*(k^2 + alpha^2*cos(phi)^2*sin(theta)^2 + alpha^2*sin(phi)^2*sin(theta)^2)^(1/2))/2 + k^2/2 + (alpha^2*cos(phi)^2*sin(theta)^2)/4 + (alpha^2*sin(phi)^2*sin(theta)^2)/4)^(3/2)*8i)/(alpha^3*cos(phi)^3*sin(theta)^3*1i - alpha^3*sin(phi)^3*sin(theta)^3 + alpha^3*cos(phi)*sin(phi)^2*sin(theta)^3*1i - alpha^3*cos(phi)^2*sin(phi)*sin(theta)^3) + (cos(theta)*(alpha^2*cos(theta)^2 + 4*k^2 + 2*alpha^2*cos(phi)^2*sin(theta)^2 + 2*alpha^2*sin(phi)^2*sin(theta)^2))/(alpha^2*cos(phi)^3*sin(theta)^3 + alpha^2*sin(phi)^3*sin(theta)^3*1i + alpha^2*cos(phi)*sin(phi)^2*sin(theta)^3 + alpha^2*cos(phi)^2*sin(phi)*sin(theta)^3*1i) + (2*(alpha^2*cos(theta)^2 + 4*k^2 + 2*alpha^2*cos(phi)^2*sin(theta)^2 + 2*alpha^2*sin(phi)^2*sin(theta)^2)*((alpha^2*cos(theta)^2)/4 - (k*(k^2 + alpha^2*cos(phi)^2*sin(theta)^2 + alpha^2*sin(phi)^2*sin(theta)^2)^(1/2))/2 + k^2/2 + (alpha^2*cos(phi)^2*sin(theta)^2)/4 + (alpha^2*sin(phi)^2*sin(theta)^2)/4)^(1/2))/(alpha^3*cos(phi)^3*sin(theta)^3 + alpha^3*sin(phi)^3*sin(theta)^3*1i + alpha^3*cos(phi)*sin(phi)^2*sin(theta)^3 + alpha^3*cos(phi)^2*sin(phi)*sin(theta)^3*1i) - (cos(theta)*((alpha^2*cos(theta)^2)/4 - (k*(k^2 + alpha^2*cos(phi)^2*sin(theta)^2 + alpha^2*sin(phi)^2*sin(theta)^2)^(1/2))/2 + k^2/2 + (alpha^2*cos(phi)^2*sin(theta)^2)/4 + (alpha^2*sin(phi)^2*sin(theta)^2)/4)*4i)/(alpha^2*cos(phi)^3*sin(theta)^3*1i - alpha^2*sin(phi)^3*sin(theta)^3 + alpha^2*cos(phi)*sin(phi)^2*sin(theta)^3*1i - alpha^2*cos(phi)^2*sin(phi)*sin(theta)^3);
c1=-((alpha^2*cos(theta)^2)/2 + 2*k^2 + (alpha^2*cos(phi)^2*sin(theta)^2)/2 + (alpha^2*sin(phi)^2*sin(theta)^2)/2)/(2*((alpha*k*cos(phi)*sin(theta))/2 - (alpha*k*sin(phi)*sin(theta)*1i)/2)) + (2*((alpha^2*cos(theta)^2)/4 - (k*(k^2 + alpha^2*cos(phi)^2*sin(theta)^2 + alpha^2*sin(phi)^2*sin(theta)^2)^(1/2))/2 + k^2/2 + (alpha^2*cos(phi)^2*sin(theta)^2)/4 + (alpha^2*sin(phi)^2*sin(theta)^2)/4))/(alpha*k*cos(phi)*sin(theta) - alpha*k*sin(phi)*sin(theta)*1i);
d1=1;

m=sqrt(a1*conj(a1)+b1*conj(b1)+c1*conj(c1)+d1*conj(d1));

u=1/m*[a1;b1;c1;d1];
rho=u*ctranspose(u);
y=kron(sigma2,sigma2); 
rhofinal=rho*y*ctranspose(rho)*y;

q=eig(rhofinal)

Moreover, is there any way to call this eigen vector for next operation without expressing its large form. Here a1,b1,c1 and d1 are basically components of a particular eigen vectors of matrix h. I pasted it from command plate

[Walter Roberson]

If I respond, are you going to delete my contributions later, like you did in the 3d plotting question?

[AVM]

@walter: Thanks for replying me.

Apologize for my nonsense activities.. I am really new user this matlab. I thought that as I asked you people so many times regarding my corrections and rectifications, so better to clear my all these garbage kind of thing so that some important and summary conversation remain. But I really don't know that it actually causes some harmful to you.

Pl help me first to recover all that deleted comment by me. I am really feeling a guilty.

[AVM]

@walter: PL help me to recover that conversation first if possible. Otherwise, I always feel ashame to ask you further.

[Walter Roberson]

It cannot be recovered.

[AVM]

@walter : I am really sorry. This happens beyond my knowldege. I am promising you that next time such kind of nonsense would not happen.

I must be carefull about all these options. Pl pardon me if possible.

[AVM]

@walter : Pl help me to solve the issue.

I would like to get eigenvalues of matrix in my code given below, the matlab taking so much time as if never ending processws. Pl somebody hepl me how to find the eigenvalues quickly.

clc
close all
syms 'theta' 'alpha' 'phi' k

%% sigma matrices are here

sigma1=[0 1;1 0];
sigma2=[0 -1i;1i 0];
sigma3=[1 0;0 -1];
sigmap=1/2*(sigma1+1i*sigma2);
sigmam=1/2*(sigma1-1i*sigma2);

%% unit vetor on sphere
n=[sin(theta)*cos(phi);sin(theta)*sin(phi);cos(theta)];

identity1=eye(2);
identity2=identity1;

p1=sigma1*sin(theta)*cos(phi)+sigma2*sin(theta)*sin(phi)+sigma3*cos(theta);
p2=kron(sigmap,sigmap);
p3=kron(sigmam,sigmam);

h=1/2*alpha*kron(p1,identity2)+k*(p2+p3);

[V,L]=eig(h);

%% The componens of one of the eigen vectors of h above pasted from command palte

a1=(4*((alpha^2*cos(theta)^2)/4 - (k*(k^2 + alpha^2*cos(phi)^2*sin(theta)^2 + alpha^2*sin(phi)^2*sin(theta)^2)^(1/2))/2 + k^2/2 + (alpha^2*cos(phi)^2*sin(theta)^2)/4 + (alpha^2*sin(phi)^2*sin(theta)^2)/4)^(3/2))/(alpha^2*k*cos(phi)^2*sin(theta)^2 + alpha^2*k*sin(phi)^2*sin(theta)^2) - (4*k^2*cos(theta) + alpha^2*cos(theta)^3 + alpha^2*cos(phi)^2*cos(theta)*sin(theta)^2 + alpha^2*cos(theta)*sin(phi)^2*sin(theta)^2)/(2*(alpha*k*cos(phi)^2*sin(theta)^2 + alpha*k*sin(phi)^2*sin(theta)^2)) + (2*cos(theta)*((alpha^2*cos(theta)^2)/4 - (k*(k^2 + alpha^2*cos(phi)^2*sin(theta)^2 + alpha^2*sin(phi)^2*sin(theta)^2)^(1/2))/2 + k^2/2 + (alpha^2*cos(phi)^2*sin(theta)^2)/4 + (alpha^2*sin(phi)^2*sin(theta)^2)/4))/(alpha*k*cos(phi)^2*sin(theta)^2 + alpha*k*sin(phi)^2*sin(theta)^2) - ((alpha^2*cos(theta)^2 + 4*k^2 + alpha^2*cos(phi)^2*sin(theta)^2 + alpha^2*sin(phi)^2*sin(theta)^2)*((alpha^2*cos(theta)^2)/4 - (k*(k^2 + alpha^2*cos(phi)^2*sin(theta)^2 + alpha^2*sin(phi)^2*sin(theta)^2)^(1/2))/2 + k^2/2 + (alpha^2*cos(phi)^2*sin(theta)^2)/4 + (alpha^2*sin(phi)^2*sin(theta)^2)/4)^(1/2))/(alpha^2*k*cos(phi)^2*sin(theta)^2 + alpha^2*k*sin(phi)^2*sin(theta)^2);
b1=-(((alpha^2*cos(theta)^2)/4 - (k*(k^2 + alpha^2*cos(phi)^2*sin(theta)^2 + alpha^2*sin(phi)^2*sin(theta)^2)^(1/2))/2 + k^2/2 + (alpha^2*cos(phi)^2*sin(theta)^2)/4 + (alpha^2*sin(phi)^2*sin(theta)^2)/4)^(3/2)*8i)/(alpha^3*cos(phi)^3*sin(theta)^3*1i - alpha^3*sin(phi)^3*sin(theta)^3 + alpha^3*cos(phi)*sin(phi)^2*sin(theta)^3*1i - alpha^3*cos(phi)^2*sin(phi)*sin(theta)^3) + (cos(theta)*(alpha^2*cos(theta)^2 + 4*k^2 + 2*alpha^2*cos(phi)^2*sin(theta)^2 + 2*alpha^2*sin(phi)^2*sin(theta)^2))/(alpha^2*cos(phi)^3*sin(theta)^3 + alpha^2*sin(phi)^3*sin(theta)^3*1i + alpha^2*cos(phi)*sin(phi)^2*sin(theta)^3 + alpha^2*cos(phi)^2*sin(phi)*sin(theta)^3*1i) + (2*(alpha^2*cos(theta)^2 + 4*k^2 + 2*alpha^2*cos(phi)^2*sin(theta)^2 + 2*alpha^2*sin(phi)^2*sin(theta)^2)*((alpha^2*cos(theta)^2)/4 - (k*(k^2 + alpha^2*cos(phi)^2*sin(theta)^2 + alpha^2*sin(phi)^2*sin(theta)^2)^(1/2))/2 + k^2/2 + (alpha^2*cos(phi)^2*sin(theta)^2)/4 + (alpha^2*sin(phi)^2*sin(theta)^2)/4)^(1/2))/(alpha^3*cos(phi)^3*sin(theta)^3 + alpha^3*sin(phi)^3*sin(theta)^3*1i + alpha^3*cos(phi)*sin(phi)^2*sin(theta)^3 + alpha^3*cos(phi)^2*sin(phi)*sin(theta)^3*1i) - (cos(theta)*((alpha^2*cos(theta)^2)/4 - (k*(k^2 + alpha^2*cos(phi)^2*sin(theta)^2 + alpha^2*sin(phi)^2*sin(theta)^2)^(1/2))/2 + k^2/2 + (alpha^2*cos(phi)^2*sin(theta)^2)/4 + (alpha^2*sin(phi)^2*sin(theta)^2)/4)*4i)/(alpha^2*cos(phi)^3*sin(theta)^3*1i - alpha^2*sin(phi)^3*sin(theta)^3 + alpha^2*cos(phi)*sin(phi)^2*sin(theta)^3*1i - alpha^2*cos(phi)^2*sin(phi)*sin(theta)^3);
c1=-((alpha^2*cos(theta)^2)/2 + 2*k^2 + (alpha^2*cos(phi)^2*sin(theta)^2)/2 + (alpha^2*sin(phi)^2*sin(theta)^2)/2)/(2*((alpha*k*cos(phi)*sin(theta))/2 - (alpha*k*sin(phi)*sin(theta)*1i)/2)) + (2*((alpha^2*cos(theta)^2)/4 - (k*(k^2 + alpha^2*cos(phi)^2*sin(theta)^2 + alpha^2*sin(phi)^2*sin(theta)^2)^(1/2))/2 + k^2/2 + (alpha^2*cos(phi)^2*sin(theta)^2)/4 + (alpha^2*sin(phi)^2*sin(theta)^2)/4))/(alpha*k*cos(phi)*sin(theta) - alpha*k*sin(phi)*sin(theta)*1i);
d1=1;

m=sqrt(a1*conj(a1)+b1*conj(b1)+c1*conj(c1)+d1*conj(d1));

u=1/m*[a1;b1;c1;d1];
rho=u*ctranspose(u);
y=kron(sigma2,sigma2);
rhofinal=rho*y*ctranspose(rho)*y;

q=eig(rhofinal)

Second thing , is there any way to call the partulcar eigen vector of a given matirx. Because, here by the expression of a1,b1,c1 and d1 I actually have written the components of one of the eigenvector of h from command plate. But I think it is not an efficient way to write.

Pl help to do that.

[Walter Roberson]

"is there any way to call the partulcar eigen vector of a given matirx."

That question is not meaningful as phrased.

Use the two output form of eig() and the eigenvectors will be the columns of the first output. You can index them directly by column number.

[AVM]

@walter: Thanks for your reply. Sorry, I didn't get your point. I know to calculate eigenvalues and eigvenvector by using [V,L]=eig(). Pl show me that in code form.

[Walter Roberson]

[V,L] = eig(h);
V(:,3)    %extract the third eigenvector

[AVM]

Thanks a lot..

[AVM]

@Walter: Here I try to calculate the concurrenec for qubit system.

Concurrence ©= max{0, e1-e2-e3-e4} where e1,e2,e3,e4 are the eigenvalues of the density matrix of R, the form given in my code.

Here e1>=e2>=e3>=e4.

Pl help me to get the expression for concurrence.

clc
close all
syms 'theta' 'alpha' 'phi' k


%% sigma matrices are here
sigma1=[0 1;1 0];
sigma2=[0 -1i;1i 0];
sigma3=[1 0;0 -1];
sigmap=1/2*(sigma1+1i*sigma2);
sigmam=1/2*(sigma1-1i*sigma2);
%% unit vetor on sphere
n=[sin(theta)*cos(phi);sin(theta)*sin(phi);cos(theta)];

identity1=eye(2);
identity2=identity1;

p1=sigma1*sin(theta)*cos(phi)+sigma2*sin(theta)*sin(phi)+sigma3*cos(theta);
p2=kron(sigmap,sigmap);
p3=kron(sigmam,sigmam);

h=1/2*alpha*kron(p1,identity2)+k*(p2+p3);

[V,L]=eig(h);

%%state vector
x1=V(:,1);

m=sqrt(dot(x1,x1));

%% Normalized state vector

x2=x1/m;

%% Density matrix
rho=x2*ctranspose(x2);
y=kron(sigma2,sigma2); 


function c=concurrence(rho);      
    R=rho*y*ctranspose(rho)*y;
  % Real is needed since MATLAB always adds a small imaginary part ...
  e=real(sqrt(eig(R)));
  e=-sort(-e);
  c=max(e(1)-e(2)-e(3)-e(4),0)
end

 
  

I need the expression of c, but here nothing is coming out.

[AVM]

@walter: How to resume the Matlab programming when it is pause situation?? And how to stop it when it is showing busy below without closing Matlab software itself?

[Walter Roberson]

What is a "pause situation" ?

If you are stopped in the debugger and wish to continue then

dbcont

If it is busy and you have the editor open, in more recent versions you can click on the Pause button . It might take some time before you get the command line control back.

If the delay for Pause is getting too long, you can control-C in the command window. There are some circumstances under which it can take a fair time to get the command line control back anyhow.

If you managed to use up all system memory and you are swapping to disk, then most often you should use your operating system task killer to kill MATLAB .

[AVM]

Thanks...Are the latest version of matlab little bit faster than 2018a? I am using 2018a. It is taking so much time for even 'simplify(f)' kind of operation and also for any complicated mathematical work.

[Walter Roberson]

There is no faster version of symbolic simplification, except possibly some fairly small gains in performance in newer versions.

Symbolic work is done using compiled libraries that are not written in MATLAB itself, so improvements in the MATLAB Execution Engine that have been made do not improve the symbolic engine.

If you are doing numeric integration on a system that has only one free variable, then use vpaintegral() instead of int() -- but integral() will likely be faster than vpaintegral()

[AVM]

Thanks for your information.

Answer 2

[AVM]

@Walter: Thanks. I am trying that

[AVM]

@Strider: Thanks for your reply.

Answer 3

[AVM]

@walter: How do I can normalize an eigen vector in Matlab? PL tell me what is the corresponding command for the normalization?

[Walter Roberson]

I do not know.

[Walter Roberson]

V_normalized = V ./ sqrt(sum(V.^2));               %requires R2016b or newer

[AVM]

Thanks a lot again.

[AVM]

@walter: With this normalized command when I try to run the following progamming, it is taking so much time as near to 1 hour and the pc start to hang frequently. I had forcefully quit that. You, pl have a look on my code and suugest me what to do. I am using matlab R2018a version.

clc;clear;

syms theta  phi  g 

%% h is 4*4 matrix

h=[ cos(theta) 0 sin(theta)*exp(-1i*phi) g
    0 cos(theta) 0 sin(theta)*exp(-1i*phi)
    sin(theta)*exp(1i*phi) 0 -cos(theta) 0
    g sin(theta)*exp(1i*phi) 0 -cos(theta) ];

%% a1,b1,c1 and d1 are the components of first eigenvector of h; these are pasted from command plate

%a1=(cos(theta)^2 + sin(theta)^2 + g^2/2 - (g*(4*sin(theta)^2 + g^2)^(1/2))/2)^(3/2)/(g*sin(theta)^2) - (cos(theta)^3 + cos(theta)*sin(theta)^2 + g^2*cos(theta))/(g*sin(theta)^2) + (cos(theta)*(cos(theta)^2 + sin(theta)^2 + g^2/2 - (g*(4*sin(theta)^2 + g^2)^(1/2))/2))/(g*sin(theta)^2) - ((g^2 + cos(theta)^2 + sin(theta)^2)*(cos(theta)^2 + sin(theta)^2 + g^2/2 - (g*(4*sin(theta)^2 + g^2)^(1/2))/2)^(1/2))/(g*sin(theta)^2);
%b1=-(exp(-phi*1i)*(cos(theta)^2 + sin(theta)^2 + g^2/2 - (g*(4*sin(theta)^2 + g^2)^(1/2))/2)^(3/2))/sin(theta)^3 + (exp(-phi*1i)*(cos(theta)^2 + 2*sin(theta)^2 + g^2)*(cos(theta)^2 + sin(theta)^2 + g^2/2 - (g*(4*sin(theta)^2 + g^2)^(1/2))/2)^(1/2))/sin(theta)^3 - (exp(-phi*1i)*cos(theta)*(cos(theta)^2 + sin(theta)^2 + g^2/2 - (g*(4*sin(theta)^2 + g^2)^(1/2))/2))/sin(theta)^3 + (exp(-phi*1i)*cos(theta)*(cos(theta)^2 + 2*sin(theta)^2 + g^2))/sin(theta)^3;
%c1= -(g^2*exp(phi*1i) + exp(phi*1i)*cos(theta)^2 + exp(phi*1i)*sin(theta)^2)/(g*sin(theta)) + (exp(phi*1i)*(cos(theta)^2 + sin(theta)^2 + g^2/2 - (g*(4*sin(theta)^2 + g^2)^(1/2))/2))/(g*sin(theta));
%d1=1;
%m=sqrt(a1*conj(a1)+b1*conj(b1)+c1*conj(c1)+d1*conj(d1));

%% normalized first eigen vector of h
%u=1/m*[a1;b1;c1;d1]; 



%% differently 

[V,L]=eig(h);

u=V(:,1)./sqrt(sum(V(:,1).^2)); %% this call aslo for normalised eigenvector of h

x=diff(u,phi);


r=dot(u,x);
assume(theta>=0);
assume(phi>=0);
r=simplify(r,'Steps',100);

f(theta,g)=1i/pi*int(r,phi,0,2*pi);
f=simplify(f,'Steps',100);

ffcn = matlabFunction(f);

theta = linspace(0.001,4, 30);
g = linspace(.001,2, 30);
[Th,G] = ndgrid(theta, g);
%F=abs(ffcn(Th,Al));
F=ffcn(Th,G);

%max(imag(F(:)))
%min(imag(F(:)))
figure
meshc(Th,G, F)
colormap(cool)
grid on
xlabel('\bf\theta','FontSize',14)
ylabel('\bf\alpha','FontSize',14)
zlabel('\bf\itf','FontSize',14)












    

But when I normalize the state just by activiting a1,b1, c1 and d1, then the graph is appearing withhin few min..But I wnated to avoid those large expression of a1,b1, c1 and d1. Pl suggest me what to do.

Thanks.