You are plotting a function of two variables, so you need to use a plot that accommodates them, and that your function is symbolic.
Try this:
syms a e
A = [besselj(0,-1i*a*e), bessely(0,-1i*a*e); ...
besselj(0,-1i*a), bessely(0,-1i*a)];
b = [-1; -1];
X = linsolve(A, b);
disp(X)
a = linspace(5,9);
e = linspace(0.05, 0.15);
fsurf(real(X), [5 9, 0.05 0.15])
xlabel('a')
ylabel('e')
Experiment to get the result you want.
