I want to be able to compute and graph the divergence of a three dimensional vector field. I have tried the following and the resulting graph is as attached.
clc
syms x y z
f = input("Enter the 3D vector function as [f1, f2, f3]");
div(x, y, z) = divergence(f, [x, y, z]);
p(x, y, z) = f(1); q(x, y, z) = f(2); r(x, y, z) = f(3);
x = linspace(-4, 4, 20); y = x; z = x;
[X, Y, Z] = meshgrid(x, y, z);
U = p(X, Y, Z); V = q(X, Y, Z); W = r(X, Y, Z);
hold on
quiver3(X, Y, Z, U, V, W, 1.5);
axis on
hold off
title("Vector Field of F(x, y, z) = [f1, f2, f3]");
and the figure output for $(x, 2y^2, 3z^3)$ -
Now obviously I would like for the figure to be represented in 3D, and I did use quiver3 for the purpose.
Surely there's something iffy with my code. An elaborate explanation about my mistakes here would be much appreciated.
I must state that I copied the code from what was taught to us for computing the divergence of 2D functions and simply went about making changes to accommodate the third variable. Also, I'm a total noob at MATLAB. My professor doesn't help my situation. Their forums have been helpful in understanding.

Your function is $f({\bf x}) = (x, 2 y^2, 3 z^3)$ so $\nabla f({\bf x}) = (1, 4 y, 9 z^2)$.
Here is a plot of the divergence (as a scalar color) atop the vector function (rendered in Mathematica):