Use the Symbolic Toolbox
a=-0.5*1i;
b=-0.8*1i;
c=0.3*1i;
[X,T]=meshgrid(-20:.5:20,-20:.5:20);
syms x t
D=2*1i*a.*x-1i*t./a;
E=2*1i*b.*x-1i*t./b;
F=2*1i*c.*x-1i*t./c;
X1=exp(D)+1i.*(exp(-D));
X2=exp(E)+1i.*(exp(-E));
X3=exp(F)+1i.*(exp(-F));
Y1=exp(D)-1i.*(exp(-D));
Y2=exp(E)-1i.*(exp(-E));
Y3=exp(F)-1i.*(exp(-F));
A=b.*X1.*Y2.*(a^3.*X3-c^3.*Y3)+X2.*(-a*b^3.*X3.*Y1+c.*((-a^3+b^3).*X1+a*c^3.^Y1).*Y3);
B=a*b*(-a+b).*X3.*Y1.*Y2+c.*(a.*(a-c).*X2.*Y1+b.*(-b+c).*X1.*Y2).*Y3;
C=(-a^3+b^3)*c.*X3.*Y1.*Y2+(b.*(a^3-c^3).*X2.*Y1+a*(-b^3+c^3).*X1.*Y2).*Y3;
D=a.*(a-c)*c.*X1.*X3.*Y2+b.*X2.*(c*(-b+c).*X3.*Y1+a*(-a+b).*X1.*Y3);
q=((1i/2).*((B.*diff(A,x)-A.*diff(B,x))./(B.*B)))+((1i/2).*((D.*diff(C,x)-C.*diff(D,x))./(D.*D)));
fprintf('optimization takes about 2 minutes, please wait\n');
tic;
qfun = matlabFunction(q, 'vars', {x, t}, 'file', 'qfun.m', 'Optimize', true);
toc
fprintf('done optimization\n');
tic
qn = qfun(X, T);
toc
subplot(2,2,1)
surf(X, T, real(qn), 'edgecolor', 'none');
title('real part')
subplot(2,2,2)
surf(X, T, imag(qn), 'edgecolor', 'none');
title('imaginary part')
subplot(2,2,[3 4])
surf(X, T, abs(qn), 'edgecolor', 'none');
title('absolute value')
A fair number of the entries come out as NaN. Some of those would not come out as NaN if you processed completely in the symbolic toolbox rather than using matlabFunction,
%this step takes a couple of minutes
qn = double( vpa(subs(q, {x,t}, {X, T})) );
