It may surprise you that only a sphere can be treated as a point mass to calculate the gravitational field near the surface applying Newton's law of universal gravitation. Other shapes can only be treated as point sources far away from them.
Calculating the gravitational field of a cylinder involves complicated integrals for the general case, that I am not able to solve analytically. I have found analytical solutions for special cases, e.g. the field along the axis of the cylinder . I could not find an analytical solution for any point around the cylinder. So I decided to implement a numerical solution.
The following graphic shows the Numerical calculation of the gravitational acceleration or gravitational potential around a cylinder. It is computed by decomposing the cylinder into many small elements and calculating for each field cell location the gravitatinoal contribution of all the elements of the cylinder. Each individual cylinder element is treated as a point mass, where the mass and center of mass depends on the individual shape of an element, see Center of Mass of a Ring Section Element.
The gravitational field near the cylinder surface does not point to the center of mass of the cylinder (black cross). In fact there is no such thing as the center of gravity for a cylinder.
Very near the surface the gravity vectors are almost vertical to the surface for some distance from the axis and the acceleration/potential are not changing much. In fact, for an infinite plane, the gravity vectors will point exactly perpendicular to the surface and the magnitude is the same everywhere, see Gravity on an infinite Flat Earth Plane for a derivation.
The analytically calculated acceleration along the axis of the cylinder is in good concordance with the numerically calculated value. Only very near the surface and inside the cylinder are some differences, depending on the choosen size of the cylinder elements. The smaller the element size, the better the concordance.
Width increasing distance from the cylinder, the gravity field approaches the field of a point mass with the center of gravity at the center of mass of the cylinder and the same mass as the cylinder. So in this case the cylinder can be regarded as a sphere or point mass of the same mass.
For a distance from the surface of about the element size or more the numerical approximation is very accurate. This can be confirmed by changing the element size and calculating the field for the same point.
The cylinder is decomposed into many small elements, see Cylinder Decomposition. For the computed field, a cross section plane is choosen and subdivided into a grid of cells. Because the cylinder is axial symmetric, this plane shows the field values for all points around the cylinder.
The volume, mass and center of mass of each cylinder element is computed according to the shape of the element. Then for each grid cell the gravitational acceleration and gravitational potential due to each cylinder element is calculated according to Newton's law of gravitation. So in this calculations each cylinder element is approximated by a point mass of the mass of the cylinder element.
The contributions of all cylinder elements are summed over by vector addition to get the final field value and direction of acceleration for one grid cell. The process is repeated for each grid cell of the cross section plane. The grid cell is colored according to the field strength and the direction of acceleration is indicated by a short line.
To cross check the validity of the computation, the analytical solution along the cylinder axis is computed as well and displayed. Additionally the field for a point mass or sphere of the same mass as the cylinder is computed.
Depending on the Radius and Height the cylinder is decomposed into small elements with an approximate size as entered in Element Size. The following graph shows how this decomposition is executed. Each element is calculated such, that the height, the radial size and the outer border size of each element is between Element Size and 2 times Element Size.
To calculate the gravitational influence of a single element, we have to know its center of mass. To calculate the center of mass of an element, we have to calculate the plane that divides the ring section element into 2 equal areas. The center of mass is then the intersection of this plane and the other perpendicular symmetry planes.
The blue and yellow areas are equal in size, see Calculation of Center of Mass Line. The red line mark shows the position of the center of mass for a triangle with 2 sides r. If r = dr and the angle φ is small, the divider line and the mark coinside.
The center of mass of any cylinder element in the radial direction of the cylinder is calculated with the following equation.
This equation is only valid for a certain maximal value of
The blue and yellow area have to be the same. We have to figure out at which distance x from the center of the circle this is the case. The blue area A is the pie area PA with radius r minus the triangle area T with height x and angle φ. The yellow area B is the triangle area T minus the pie area PB with radius r − Δr. We need to find the value x for which the areas A and B are equal:
The area of the triangle T depends on its height x and half the base length s, which depends on the height x and the angle φ:
The areas of the pie sections of the circle are:
Putting this together we get:
Now we can solve for x:
There is a relative simple analytical equation for the gravitational acceleration of a cylinder along its symmetry axes. The result of the analytical equation is displayed in the App at the lower right corner. It should be close to the value displayed at the cursor position, if the cursor is on the vertical y-axis. Near the surface and inside the cylinder the value from the numerical calculation and the analytical value can differ slightly due to the limited resolution of the numerical solution.
Analytical solutions also exist for the potential of a cylinder but are not included here.
The gravitational field of a cylinder along its axis is: 
Depending on whether the point on the y-axis is outside or inside the cylinder we have to provide different bounds b for the integral. The y values are measured from the surface and always positive.
If we do the integral for both cases we get: