contour plot of charge density over electric field

close all;
clear all;
syms x y ep0 %defined x,y,ep0 as symbols
k=10^-9/(36*pi); %value of epsilon
E=[x^4 y^3]; %expression for electric field
rho=ep0*divergence(E,[x,y]); %expression for rho
rho=subs(rho,ep0,k);
v=-2.5:0.1:2.5;
[xP,yP]=meshgrid(v);
Ex=subs(E(1),x,xP);
Ey=subs(E(2),y,yP);
rhoPlot=double(subs(rho,{x,y},{xP,yP}));
figure(1)
quiver(xP,yP,Ex,Ey);
hold on
contour(xP,yP,rhoPlot)
title('contour plot of charge density over electric field')
xlabel('x');
ylabel('y');