It is not easy because there is no closed form expression for the median and quartiles. For the median you can find [several approximate formula's](https://en.wikipedia.org/wiki/Gamma_distribution#Median_approximations_and_bounds). For the quartiles you could use R's approximation of the quantile function, `qgamma`, and then use an iterative algorithm that finds $k$ and $\theta$ that suit your target.

One way to make such algorithm can be by optimizing a cost function




    f = function(par, targetmedian, targetIQR) {
      k = par[1]
      theta = par[2]
    
      median = theta*k*(1-1/9/k)^3
      ### alternatively compute the median like below
      #median = qgamma(0.5,shape = k,scale = theta)

      IQR = qgamma(0.75,shape = k,scale = theta) - qgamma(0.25,shape = k,scale = theta)

      error = (median-targetmedian)^2+(IQR-targetIQR)^2
      return(error)
    }
    
    ### optim function below gives result k = 2.905083 theta = 1.160847
    p2 = optim(par  = c(1,1), f, targetmedian = 3, targetIQR = 2.5)$par

    ### check median and IQR
    qgamma(0.5, shape = p2[1], scale = p2[2]) ### 2.994302
    qgamma(0.75,shape = p2[1], scale = p2[2])-qgamma(0.25,shape = p1[1],scale = p1[2]) # 2.494646