variance function on a mean polynomial (nlme R package)-CodePudding

I want to fit a variance function that involves polynomial terms of the mean

$https://latex.codecogs.com/svg.image?\lambda(\mu_i)=exp(\alpha_1\mu_i \alpha_2\mu_i^{2})$

I know that I can fit the variance function

$https://latex.codecogs.com/svg.image?\lambda(\mu_i)=exp(\alpha_1\mu_i)$

using varExp() (without arguments). Reading the documentation of varExp() it defaults to ~ fitted(.), however using varExp(~ fitted(.)) gives an error, therefore varExp(~ fitted(.)^2) is not working.

Here is a code example in which I am working:

library(pacman)
p_load(tidyverse)
p_load(MASS)
p_load(nlme)
set.seed(123)
n<-1000

z0<-data.frame(days=runif(n,-3,3))

z0_mat <- model.matrix(~ days, z0)
betas<-c(10,.85)
a0<-1.5
b0<-8
b1<- -0.2
x_target<-60
alpha<-(log(b0)-log(a0))/x_target
z0$weight <- z0_mat%*