It is unclear in your question what "manual calculation" excludes, and your comment that you cannot use "advanced functions" is also not very helpful. In any case, fitting an ARMA model via maximum-likelihood estimation is an optimisation problem where you need to maximise a function over a set of parameters. For example, the log-likelihood function in a stationary Gaussian ARMA is:
$$\ell_{x}(\mu,\boldsymbol{\phi},\boldsymbol{\theta}) = - \frac{1}{2} \ln | \boldsymbol{\Sigma}(\boldsymbol{\phi},\boldsymbol{\theta})| + (\mathbf{x} - \mu \boldsymbol{1})^\text{T} \boldsymbol{\Sigma}(\boldsymbol{\phi},\boldsymbol{\theta})^{-1} (\mathbf{x} - \mu \boldsymbol{1}),$$
where the covariance matrix $\boldsymbol{\Sigma}(\boldsymbol{\phi},\boldsymbol{\theta})$ depends on the parameters $\boldsymbol{\phi}$ and $\boldsymbol{\theta}$ according to the auto-covariance function for the ARMA model. This function has two terms. The second is a standard sum-of-squares term, but the first is a more complicated term involving the logarithm of the determinant of the covariance matrix.
The exact MLE method has critical point equations that cannot be put into closed form, so this would entail the use of iterative methods (e.g., Newton-Raphson iteration) to find the maximising values. If you are willing to deviate slightly from the exact MLE and use the partial likelihood function ---excluding the logarithmic term--- this gives MLEs that can be obtained as standard WLS estimates. Once you estimate the parameters in the model you can make forecasts as point-estimates by substituting the parameter estimates.
It is certainly possible to program this optimisation problem "manually", in the sense that you can directly program an iterative procedure to optimise the above function, without using pre-programmed optimisation procedures. It would be quite laborious, but it could probably done in a few hours.
