How to plot contour as a function of radius and height?

I have written a code for magnetic field on surface as a function of radius and height.. I want graph of variation in magnetic field on the surface as the magnetic lines appear i.e., diverging.. The code is here: %magnetic field normal to the surface

clc;
clear all;
close all;

hob=2.5*10^3; 

h=200000:100000:5000000;

c=3*10^8; %m/s velocity

len1=length(h);

bsquare=10^10:(10^15-10^10)/(len1-1):10^15; %m2

b=bsquare.^0.5;

len2=length(b);

area=bsquare.*pi;

ho=zeros(len1,len2); %mag field

fh=zeros(len1,len2); %electron gyro freq

wh=zeros(len1,len2); %ang freq

m=9.1094*10^-28; %gram electron mass

e=1.6*10^-19; %charge

for i=1:len1

      for j=1:len2

      ho(i)=hob*(1-(h(i)/((h(i)^2)+bsquare(j))^0.5));

      wh(i)=e*ho(i)/(m*c);

      fh(i)=wh(i)/(2*pi);

      end

end

figure;

t=[h:b];

contour3(b,h,ho,128);%logspace(-2, 3, 20));

% plot3(b,h,ho);

xlabel('x-radius');ylabel('y-height');zlabel('z-magnetic field');

grid on;

hold on;

[Saumitra Lalsare]

I got the answer very well...

The code goes as follows:

clc;
clear all;
close all;

hob=2.5*10^3; %oe

%h=input('h= ');

%b=input('b= ');

h=20000000:10000000:100000000;  % height

b=10^6:10^7:10^9;               % radius

n=length(h);m=length(b);

ho=zeros(m,n);                  % magnetic field

for i=1:n

      for j=1:m

          ho(j,i)=hob*(1-(h(i)/((h(i)^2+b(j)^2)^0.5)));

      end

end

figure;

subplot(1,2,1);

surf(h,b,ho);title('Magnetic field');xlabel('x-height');ylabel('y-
radius');zlabel('z-magnetic field');

subplot(1,2,2);

contour(h,b,ho,2048);title('Magnetic field');xlabel('x-height');ylabel('y-radius');%zlabel('magnetic field');